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Abstract 

Two  stand-level  hydrologic  computer  models  have  been  developed 
to  study  water  flow  through  western  coniferous  forest  ecosystems. 
The  models,  H20TRANS  and  DAYTRANS,  are  for  hourly  and  daily 
timesteps,  respectively,  and  use  routine  meteorological  data.  Re- 
quired parameters  include  leaf  area  index,  sapwood  basal  area,  and 
soil  storage  capacity.  The  models  incorporate  rates  of  snowmelt, 
precipitation,  canopy  interception,  and  litter  and  soil  evaporation. 
Primary  model  outputs  are  transpiration,  soil  moisture  depletion,  sub- 
surface outflow,  and  tree  water  stress  development  as  measured  by 
leaf  water  potential  and  leaf  conductance.  Complete  documentation 
of  all  equations  is  presented,  as  are  the  results  from  an  initial  valida- 
tion on  lodgepole  pine  at  the  Fraser  Experimental  Forest  in  Colorado. 
A  sensitivity  analysis  of  the  models  and  a  discussion  of  their  range  of 
applicability  is  also  presented. 


'Previously  with  the  Mountain  Meteorology  Project,  Rocky  Mountain  Forest  and  Range  Experi- 
ment Station,  Fort  Collins,  Colo. 
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Documentation  and  Preliminary  Validation 

of  H20TRANS  and  DAYTRANS, 
Two  Models  for  Predicting  Transpiration 
and  Water  Stress  in  Western  Coniferous  Forests 

Steven  W.  Running 


INTRODUCTION 

Water  movement  in  a  forest  stand  is  determined  by 
both  physical  and  physiological  factors.  A  forest 
ecosystem  is  diverse  in  species  composition,  physical 
site  characteristics,  and  micrometeorology.  Physical 
factors,  such  as  precipitation,  snowmelt,  and  soil 
storage  provide  the  water  supply,  while  atmospheric 
evaporative  force  generates  the  demand.  In  a  forest,  the 
tree  is  a  major  water  transport  system  from  the  soil  to 
the  atmosphere;  this  transport  system  is  controlled  by 
physiological  responses  by  the  tree  to  the  physical 
environment. 

Some  researchers  develop  stand  water  use  models 
primarily  from  a  hydrologic  perspective  (Leaf  and  Brink 
1973).  However,  recent  models  of  the  Soil-Plant- 
Atmosphere  Continuum  (SPAC)  have  blended  micro- 
meteorology,  hydrology,  and  physiology.  Recent  models 
using  this  approach  have  been  reported  by  Sinclair  et  al. 
(1976),  Tan  et  al.  (1978),  Luxmoore  et  al.  (1978),  Federer 
(1979),  and  Lohammar  et  al.  (1980). 

The  two  models  presented  in  this  paper,  H20TRANS 
and  DAYTRANS,  have  been  developed  primarily  to  esti- 
mate transpiration  and  water  stress  development  in 
western  coniferous  forests.  Earlier  development  of  these 
models  originated  primarily  from  physiological  research 
on  stomatal  control  of  conifers  under  stress  (Running  et 
al.  1975,  Waring  and  Running  1976).  More  recent  efforts 
continued  the  emphasis  on  tree-water  relations  but 
added  more  physical  processes,  such  as  snowmelt  and 
canopy  evaporation,  to  produce  a  more  accurate  model 
for  wider  application. 

The  philosophy  of  both  models  considers  a  stand  to  be 
a  collection  of  individual  trees.  This  perspective  is 
necessary  because  much  of  the  physiological  control  of 
water  movement  in  a  stand  occurs  within  the  tree.  A 
stand  is  a  defined  unit,  whereas  a  tree  is  a  biological 
unit.  Consequently,  while  some  parts  of  the  models  are 
stand  oriented,  such  as  leaf  area  index  and  canopy  light 
attenuation,  the  models  can  also  be  programmed  to 
simulate  individual  trees.  This  is  not  without  some  dif- 
ficulty, such  as  having  to  estimate  lateral  rooting 
distance  of  an  individual  tree  instead  of  defining  ground 
surface  area  as  total  stand  area. 

For  simplicity  the  models  do  not  layer  the  canopy  into 
different  compartments  or  consider  vertical  or  horizon- 
tal micrometeorological  profiles.  A  single  canopy 
average  is  used  for  meteorological  inputs  and  for  com- 


putation of  physiological  responses  (i.e.,  canopy  average 
leaf  conductance  and  leaf  water  potential).  Because  the 
models  do  not  incorporate  any  terrestrial  subunit 
stratifications,  stands  must  be  relatively  homogeneous 
topographically  for  optimum  model  performance. 

Many  of  the  functions  used  in  the  models  were 
developed  from  research  at  the  Fraser  Experimental 
Forest  in  central  Colorado,  working  on  Pinus  contorta 
Dougl.  var.  latifolia  Engelm.  (lodgepole  pine),  and  from 
research  on  Pseudotsuga  menziesii  (Mirb.)  Franco 
(Douglas-fir)  in  western  Oregon. 


MODEL  DOCUMENTATION 

H20TRANS  and  DAYTRANS  are  two  process-level 
compartment  simulation  models  that  use  climatic, 
edaphic,  and  stand  data  to  estimate  transpiration,  soil 
moisture  depletion,  and  water  stress  development  in  a 
forest  stand.  The  primary  difference  between  the  two 
models  is  that  H20TRANS  runs  on  hourly  input  climatic 
data  while  DAYTRANS  only  requires  daily  average  in- 
puts of  climatic  data.  The  most  significant  outputs  from 
the  models  are  soil  moisture  depletion,  stand  transpira- 
tion, leaf  conductance,  and  leaf  water  potential.  Com- 
partment diagrams  with  flow  linkages  are  illustrated  in 
figures  1  and  2. 

Both  models  were  constructed  using  a  general  simula- 
tion processor  similar  to  the  FLEX/REFLEX  system 
developed  in  the  IBP  Coniferous  Forest  Biome  by  White 
and  Overton  (1977)  for  large-scale  ecosystem  models. 
Application  of  the  FLEX  modeling  format  also  can  be 
found  in  the  description  of  the  CONIFER  model  of  the 
U.S.  IBP  Coniferous  Forest  Biome  (CFBMG  1979,  Swart- 
zman  1979).  The  processor  assigns  all  parameters  a 
subscripted  variable  so  that  changes  in  the  model  can  be 
made  easily.  State  variables  or  compartments  are  desig- 
nated as  X(i).  The  change  in  X(i)  from  (t)  to  (t  +  1)  is  I(i). 
Flow  between  state  variables  is  designated  F(i,  j),  where 
material  flows  from  X(i)  to  X(j).  Driving  variables,  which 
typically  input  climatic  variables  in  ecological  models, 
are  Z(i).  Constant  parameters  used  in  the  model  are  B(i). 
Finally,  intermediate  variables  used  to  compute  the 
functions  controlling  flow  in  the  model  are  G(i).  Other 
variable  notations  are  given  in  appendix  1. 

The  processor  is  organized  as  a  command  program 
that  sequentially  calls  a  set  of  subroutines  containing 
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Figure  1.— Compartment  diagram  and  flow  linkages  for  H20TRANS.  The  double-lined  boxes 

represent  summing  compartments. 
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Figure  2.— Compartment  diagram  and  flow  linkages  for  DAYTRANS.  The  double-lined  boxes 

represent  summing  compartments. 
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model  detail.  The  sequence  and  a  brief  description  of 
subroutines  called  follows. 


Subroutine 

HEADER 


ZCOMP 


PROCESS 


FLOW 


YCOMP 
PRINTER 

UPDATE 
ZERO 


ERROR 


Description 

Reads  in  initial  information  on  the 
number  of  X,  B,  G,  and  Z  variables, 
timestep,  and  other  details  that  are 
covered  in  appendix  2. 

Reads  in  the  Z  (driving)  variables 
from  a  data  file,  explained  in  appen- 
dix 3. 

Calculates  all  the  G  (intermediate) 
variables. 

Uses  the  G  variables  calculated  in 
the  PROCESS  subroutine  to  calcu- 
late the  F  (flow)  variables,  which 
transfer  water  between  X  (state) 
variables. 


tions  for  all  state  variables.  Examples  of  initial  condi- 
tions for  the  state  variables  are  given  in  the  Results  sec- 
tion. A  list  and  descriptions  of  state  variables  follow. 


Prepares  a  matrix 
values  for  printout. 


of  computed 


Labels  the  line  printer  page,  estab- 
lishes the  output  format,  and  prints 
out  all  requested  parameters. 

Updates  the  status  of  the  X  (state) 
variables  after  each  iteration  by 
adding  or  subtracting  water  as  cal- 
culated in  the  FLOW  subroutines. 

Sets  the  I,  Z,  F,  and  G  variables  to 
zero  after  the  last  iteration  in  a 
model  run.  The  last  printout  then 
gives  a  final  condition  of  the  X 
variables. 

Is  a  library  of  error  statements 
called  when  mistakes  are  made, 
which  identifies  the  subroutine  that 
caused  problems. 

H20TRANS 


STATE  VARIABLES— X(i) 

State  variables  or  compartments  are  designated  for 
parts  of  the  stand  hydrologic  system  that  have  a  water 
storage  function.  The  state  of  the  system  can  be  partial- 
ly determined  by  the  location  of  water  in  the  system  at 
any  time.  Although  the  state  variables  give  no  indication 
of  water  transfer  in  the  system,  the  change  in  water  con- 
tent of  the  state  variables  does.  This  change  is  produced 
by  the  flow  functions  described  that  couple  the  state 
variables. 

A  new  state  variable  can  be  added  to  the  model  by 
designating  a  new  subscripted  X  in  the  header  file 
(appendix  2).  The  header  file  then  reads  the  initial  condi- 


Unit 

X(l)  = 

snowpack  water  content 

cm3 

X(2)  = 

soil  surface  runoff 

cm3 

X(3)  = 

water  storage  in  forest  floor 

litter 

cm3 

X(4)  = 

available  root  zone  water  in 

Soil  Level  A 

cm3 

X(5)  = 

available  root  zone  water  in 

Soil  Level  B 

cm3 

X(6)  = 

subsurface  outflow 

cm3 

X(7)  = 

stem  sapwood  water 

storage 

cm3 

X(8)  = 

leaf  water  storage 

cm3 

X(9)  = 

transpiration 

cm3 

X(10)  = 

evaporation,  canopy  and  lit- 

ter surfaces 

cm3 

X(l):  The  snowpack  is  defined  in  terms  of  the  water 
available  for  melt  into  the  soil.  Because  growing  season 
water  availability  is  of  greatest  importance  to  this 
model,  X(l)  is  defined  when  the  snowpack  is  isothermal 
at  0°  C  or  "ripe"  just  prior  to  spring  snowmelt.  Any  ad- 
ditional snowfall  is  treated  as  general  precipitation  by 
the  model  and  is  not  added  to  X(l).  This  was  done  to 
avoid  having  to  deal  with  the  volume  of  information 
needed  to  calculate  an  energy  budget  for  the  snowpack 
as  developed  by  Leaf  and  Brink  (1973).  Although  this 
simplified  snowpack  definition  is  not  an  accurate 
estimator  of  snowpack  dynamics,  on  most  years  it 
should  be  adequate  to  follow  snowmelt  into  the  soil. 

X(2):  An  unlimited  summing  compartment  for  surface 
runoff  water  that  occurs  when  precipitation  or  snow- 
melt exceeds  the  soil  surface  infiltration  capacity. 

X(3):  The  water  content  of  the  forest  floor  litter. 
Although  this  is  a  small  component  of  the  stand  water 
balance,  the  water  content  of  the  forest  floor  litter  is  a 
critical  factor  in  litter  decomposition  rates,  an  impor- 
tant part  of  the  terrestrial  nutrient  cycle. 

X(4),  X(5):  The  soil  water  supply  available  to  the  roots 
is  partitioned  between  X(4)  and  X(5).  Rooting  depth  and 
lateral  extent  must  be  estimated  to  compute  the  total 
volume  of  soil  tapped  by  the  root  systems.  While  approx- 
imate estimation  of  X(4)  and  X(5)  for  a  stand  is  not  dif- 
ficult, exact  measurement  of  these  compartment  sizes 
for  an  individual  tree  is  laborious.  X(4)  is  defined  as  the 
upper  50%  of  the  root  zone  and  X(5)  as  the  lower  50%, 
although  this  distribution  can  easily  be  changed.  The 
main  purpose  for  splitting  the  root  zone  water  into  two 
compartments  was  to  provide  the  option  of  simulating 
competition  between  plants  of  different  rooting  depth. 

X(6):  An  unlimited  summing  compartment  for  water 
that  percolates  below  the  rooting  zone.  This  compart- 
ment only  receives  water  when  input  precipitation  or 
snowmelt  continues  after  X(3),  X(4),  and  X(5)  have  been 
filled. 

X(7):  The  water  available  for  transpiration  from  the 
stem  sapwood  of  the  trees. 
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X(8):  Water  stored  in  the  leaves  (needles)  available  for 
transpiration.  This  source  is  small  when  compared  to 
the  rest  of  the  tree  water  system,  but  it  directly  deter- 
mines leaf  water  potential  [\p{).  Therefore,  it  is  critical  in 
the  control  of  stomatal  response  and  leaf  conductance 
(k,). 

X(9):  Cumulative  transpiration  from  the  tree  or  stand. 

X(10):  Cumulative  evaporation  from  the  stand  canopy 
and  litter  surfaces. 

At  each  iteration,  the  current  status  of  each  state 
variable  X(i)  is  printed.  In  addition,  the  net  change  in 
each  state  variable  from  the  previous  iteration  is  also 
printed,  designated  I(i).  For  example,  X(9)  is  cumulative 
transpiration  since  the  beginning  of  the  simulation,  and 
1(9)  is  the  transpiration  that  occurred  in  the  previous 
1-hour  timestep. 

FLOW  FUNCTIONS— F(i,j) 

The  convention  of  F(i,j)  is  that  X(i)  is  the  state  variable 
from  which  water  is  withdrawn;  X(j)  is  the  compartment 
to  which  water  from  X(i)  is  added.  If  i  =  j,  this  denotes 
water  input  from  a  source  external  to  the  model  enter- 
ing X(i)  or  output  to  an  unmonitored  sink.  The  F(i,j)  nota- 
tion couples  the  correct  state  variables,  while  the  G 
variables  shown  are  the  calculated  magnitude  of  flow 
between  X(i)  and  X(j)  at  any  iteration.  The  equations 
used  to  calculate  the  G  variables  are  described  in  a 
subsequent  section.  A  list  and  descriptions  of  flow  func- 
tions follow. 


Unit 

F(3,3)  = 

G(50)  -  G(53) 

cm3 

F(l,3)  = 

G(51) 

cm3 

F(3,2)  = 

G(54) 

cm3 

F(3,4)  = 

G(55) 

cm3 

F(4,5)  = 

G(56) 

cm3 

F(5,6)  = 

G(57) 

cm3 

F(8,9)  = 

G(18) 

cm3 

F(7,8)  = 

G(62) 

cm3 

F(4,7)  = 

G(67) 

cm3 

F(5,7)  = 

G(68) 

cm3 

F(10,10) 

=  G(53)  +  G(58) 

cm3 

F(3,3):  Precipitation  into  the  model  from  the  exterior 
source  (driving  variable)  minus  canopy  interception  that 
does  not  reach  the  forest  floor.  Also  subtracts  litter  sur- 
face evaporation  from  X(3). 

F(l,3):  Snowmelt  input  to  the  forest  floor  litter  com- 
partment. 

F(3,2):  Surface  runoff  that  occurs  when  precipitation 
and/or  snowmelt  exceeds  the  soil  surface  infiltration 
capacity. 

F(3,4):  Water  that  flows  into  the  root  zone  Soil  A  com- 
partment, X(4),  after  the  forest  floor  litter  is  saturated 
(i.e.,  after  X(3)  is  at  capacity). 

F(4,5):  Water  that  enters  after  X(4)  has  reached 
capacity  is  cascaded  to  X(5). 

F(5,6):  When  X(3),  X(4),  and  X(5)  are  all  at  capacity, 
additional  water  input  to  the  system  is  routed  to  subsur- 
face outflow,  X(6),  from  X(5). 


F(8,9):  Moves  water  from  leaf  storage,  X(8),  to  the 
transpiration  sink,  X(9). 

F(7,8):  Water  withdrawn  from  internal  tree  sapwood 
storage  to  help  satisfy  the  water  demand  generated  by 
F(8,9)  and  the  deficit  created  in  X(8). 

F(4,7),  F(5,7):  When  a  water  deficit  is  created  in  X(7) 
by  transpiration,  F(8,9),  this  demand  is  met  by  root 
water  uptake  from  the  two  rooting  zone  soil  reserves, 
X(4)  and  X(5). 

F(10,10):  Canopy  and  litter  evaporation. 

Note  that  water  movement  into  the  soil  state  vari- 
ables, X(l)  through  X(6),  is  basically  input  or  donor  con- 
trolled. Water  withdrawal  from  the  system,  (transpira- 
tion) is  an  output  demand  with  donor  control. 

DRIVING  VARIABLES— Z(i) 

These  variables  are  input  from  an  exterior  data  file. 
Appendix  3  presents  this  data  file  format.  Precipitation 
and  soil  temperature  are  input  once  every  24  hours;  the 
other  meteorological  variables  are  entered  hourly.  This 
was  done  because  precipitation  is  rarely  recorded  more 
than  once  daily  by  most  installations,  and  soil  tempera- 
ture at  20  cm  changes  slowly.  The  24-hour  precipitation 
is  input  into  the  model  as  1/4  of  daily  total  at  hours  1,  2, 
3,  and  4  of  each  day.  This  limits  the  effects  of  short  but 
severe  thunderstorms.  In  particular,  surface  runoff, 
F(3,2),  rarely  occurs.  Precipitation  could  easily  be  input 
hourly,  and  I  recommend  this  modification,  if  hourly 
data  can  be  recorded.  The  other  variables  either  can  be 
recorded  once  each  hour,  typically  on  the  hour,  or  be 
full-hour  averages.  A  list  of  driving  variables  follows. 


Unit 

Z(l) 

=  Julian  date 

day 

Z(2) 

=  precipitation 

cm  day-1 

Z(3) 

=  air  temperature 

°C 

Z(4) 

=  relative  humidity 

% 

Z(5) 

=  soil  temperature  (20  cm  depth) 

°C 

Z(6) 

=  incoming  shortwave  radiation 

ly  min-1 

Z(7) 

=  hour  of  day 

h 

AUXILIARY  CONSTANTS— B(i) 

The  constant  parameters  are  values  entered  into  the 
model  for  use  in  various  computations;  they  do  not 
change  throughout  the  model  run.  These  can  include 
stand  physical  dimensions,  conversion  constants,  and 
empirical  coefficients  used  to  calculate  intermediate  (G) 
variables.  All  B  parameters  are  specified  and  entered 
into  the  program  by  the  HEADER  file,  explained  in  detail 
in  appendix  2.  Examples  and  derivations  of  values  used 
for  the  B  constants  are  covered  in  the  Results  section 
because  many  are  site-specific.  A  list  of  auxiliary  con- 
stants follows. 

Unit 

B(l)    =  snowpack  melt  coefficient        cm  °C_1 
B(2)    =  maximum  soil  surface 

infiltration  rate  cm  h-1 
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cmj 


cm  LAI 


cmJ 


cnr 


cm' 
cm2 


m 


cm  s" 


MPa 


B(3)    =  water  storage  capacity  of 
forest  floor  litter,  X(4) 

B(4)    =  maximum  available  water 
storage  in  root  zone  Soil 
A,  X(4) 

B(5)    =  maximum  available  water 
storage  in  root  zone  Soil  B, 
X(5) 

B(6)    =  canopy  interception 

coefficient 
B(7)    =  maximum  available  water 

storage  in  stem  sapwood, 

X(7) 

B(8)    =  maximum  available  leaf 

water  storage  in  X(8) 
B(9)    =  tree  or  stand  total  leaf 

area  (not  projected  area) 
B(10)  =  ground  surface  area 
B(ll)  =  average  midcrown  tree 

height 

B(12)  =  maximum  canopy  average 
leaf  conductance  (k,) 

B(13)  =  spring  minimum  predawn 
leaf  water  potential  (Bi/',) 
NOTE:  Throughout  the  model, 
water  potential  is  treated  as  a 
positive  value  for  mathematical 
simplicity. 

B(14)  =  leaf  osmotic  potential  (and 
stomatal  closure  threshold) 

B(15)  =  sapwood  basal  area 

B(16)  =  leaf  water  depletion 
coefficient 

B(17)  =  stem  and  soil  water 
exchange  coefficient 


INTERMEDIATE  VARIABLES— G(i) 

The  intermediate  variables  are  used  for  a  variety  of 
preliminary  calculations  in  the  model.  Because  each 
flow  variable  is  controlled  by  many  factors,  it  is  not 
possible  to  write  one  equation  representing  all  condi- 
tions. Consequently,  G  variables  are  used  to  provide 
more  manageable  intermediate  computations.  Often  the 
G  variables  will  represent  specific  processes  of  impor- 
tance in  the  model,  such  as  the  control  of  leaf  conduct- 
ance by  incoming  radiation.  Because  all  G  variables 
have  a  unique  subscript,  they  do  not  have  to  be  used  se- 
quentially. For  clarity  it  is  useful  to  group  G  variables 
that  are  performing  similar  computations.  A  brief  index 
of  G  variables  precedes  their  detailed  documentation. 


MPa 

cm2 


G(l)  -G(9): 

G(10)  -  G(19): 

G(50)  -  G(59) 
G(60)  -  G(69) 
G(70)  -  G(79) 


water  storage  and  water  content 
calculations 

leaf  conductance  and  transpiration 
calculations 

soil  water  input  calculations 
tree  water  uptake  calculations 
leaf  water  potential,  flow  resistance 
T/PT  ratio. 


Description  of  G  Variables 

Water  Storage  and  Water  Content 
G(l)  =  absolute  humidity  deficit  (g  cm"3) 
G(l)  =  S3(Z(3),  Z(4)) 

where 
Z(3)  =  air  temperature, 
Z(4)  =  relative  humidity. 

This  statement  sends  air  temperature  and  relative 
humidity  to  Function  S3.  The  function  calculates  ab- 
solute humidity  deficit  and  returns  it  to  G(l). 

G(2)  =  leaf  area  index,  LAI  (m2  m"2) 
G(2)  =  B(9)/B(10) 

where 

B(9)  =  total  stand  (or  tree)  leaf  area,  all  surfaces, 
B(10)  =  total  stand  (or  tree)  ground  surface  area. 

G(3)  =  available  water  fraction  in  X(3),  forest  floor  litter 
(cm3  cm-3) 

G(3)  =  X(3)/B(3) 

where 

X(3)  =  present  water  content  of  the  forest  floor  litter, 
B(3)  =  water  capacity  of  the  forest  floor  litter. 
Converting  available  water  in  compartments  to  a 
fraction  allows  easy  changing  of  compartment  sizes 
through  the  Header  file  without  reprogramming  each  G 
variable. 

G(4)  =  available  water  fraction  in  X(4),  the  upper  soil 
root  zone  (cm3  cm-3) 

G(4)  =  X(4)/B(4) 

where 

X(4)  =  present  water  content  of  the  upper  soil  root 
zone, 

B(4)  =  water  capacity  of  the  upper  soil  root  zone. 

G(5)  =  available  water  fraction  in  X(5),  the  lower  soil 
root  zone  (cm3  cm-3) 

G(5)  =  X(5)/B(5) 

where 

X(5)  =  present  water  content  of  the  lower  soil  root 
zone, 

B(5)  =  water  capacity  of  the  lower  soil  root  zone. 

G(6)  =  averaged  available  water  fraction  in  the  total 
soil  rooting  zone  (cm3  cm-3) 

G(6)  ='  (G(4)  +  G(5))/2 

where 

G(4)  =  available  water  fraction  in  the  upper  soil  root 
zone, 

G(5)  =  available  water  fraction  in  the  lower  soil  root 
zone. 

G(7)  =  available  water  fraction  in  X(7),  stem  storage 
(cm3  cm-3) 

G(7)  =  X(7)/B(7) 
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where 

X(7)  =  present  stem  sapwood  water  content, 
B(7)  =  maximum    available    stem    sapwood  water 
storage. 

Leaf  Conductance  and  Transpiration 

G(10)  =  predawn  leaf  water  potential,  B\j/}  (-MPa) 

G(I0)  =  AMAX1  (B(13),  0.2/G(6)  +  0.01  ■  B(ll),  0.0) 
where 

AMAX1  =  an  intrinsic  FORTRAN  function  that 
chooses  the  maximum  from  a  set  of 
values  in  the  parentheses, 

B(13)  =  spring  minimum  predawn  leaf  water  poten- 
tial, 

0.2     =  an  empirical  coefficient  used  to  generate  the 

curve  in  figure  3, 
G(6)    =  averaged  total  available  water  fraction, 
0.01    =  hydrostatic  gradient  constant  (MPa  m1), 
B(ll)  =  midcrown  tree  height. 

This  statement  calculates  the  predawn  leaf  water 
potential  (B^,)  as  a  function  of  available  root  zone  soil 
water  in  X(4)  and  X(5)  (fig.  3).  This  type  of  relation  be- 
tween Bi/',  and  soil  water  content  has  been  found 
numerous  times  in  various  plant  communities  (Sucoff 
1972,  Hinckley  and  Ritchie  1973,  Branson  et  al.  1976, 
Huzulak  1977).  The  function  basically  combines  the 
well-known  curvilinear  relationship  between  soil  water 
content  and  soil  water  potential  (Heth  and  Kramer  1975, 
Nnyamah  and  Black  1977),  and  the  linear  relationship 
between  soil  water  potential  and  predawn  leaf  water 
potential  (Hinckley  et  al.  1978,  Drew  and  Ferrell  1979, 
Sala  et  al.  1981). 

This  treatment  of  the  soil-plant  interaction  is  rather 
superficial  and  yet,  appears  adequate  for  the  purposes 
of  this  model.  Implicit  in  the  equation  is  the  moisture  ten- 
sion release  curve  of  a  well-drained  sandy  loam.  The 
equation  would  not  easily  incorporate  a  soil  of  radically 
stratified  physical  characteristics.  In  that  instance,  it 


ol  1  1  1  i  i  i  i  i  i  I 

0  0  0.2  0  4  0  6  0  8  1.0  x  100 

Soil  water  content  (%  of  capacity) 

Figure  3.— The  function  relating  soil  water  content  (percent  of  field 
capacity)  to  predawn  leaf  water  potential,  (an  estimate  of  soil 
water  potential). 


would  be  best  to  define  separate  soil  compartments 
(X(4),  X(5),  etc.)  and  treat  each  soil  horizon  individually. 

This  equation  also  assumes  root  distribution  adequate 
to  extract  water  from  the  entire  defined  rooting  zone  soil 
compartment.  Consequently,  both  the  definition  of  avail- 
able water  in  rooting  zone  soil  (X(4)  and  X(5))  and  a  cor- 
rect response  function  for  G(10)  (fig.  3)  is  necessary  for 
the  model  to  correctly  predict  the  magnitude  and  timing 
of  soil  moisture  depletion  and  tree  water  stress  develop- 
ment. For  accurate  results,  this  response  curve  (fig.  3) 
should  always  be  checked  in  a  new  area  with  field 'meas- 
urements of  predawn  leaf  water  potential  plotted 
against  soil  water  content. 

G(ll)  =  light-saturated  leaf  area  index  (m2  m~2) 
G(ll)  =  10.0  •  (1.0  -  EXP(-4.6-Z(6))) 

where 

10.0,  1.0,  -4.6  =  curve-fitting  coefficients, 

EXP  =  exponential, 

Z(6)    =  incoming  shortwave  radiation. 

Stomata  of  many  trees  do  not  fully  open  until  incoming 
shortwave  radiation  reaches  around  10%  of  full  sun- 
light (Davies  and  Kozlowski  1974,  Dykstra  1974,  Hinck- 
ley et  al.  1975,  Watts  et  al.  1976,  Hinckley  et  al.  1978). 
So,  G(ll)  calculates  how  much  leaf  area  index  is  above 
the  radiation  threshold  of  0.1  ly  min-1  at  the  current 
radiation  level  Z(6).  A  typical  Beer's  Law  exponential 
decrease  in  radiation  as  it  passes  vertically  through  the 
canopy  is  assumed  (fig.  4).  The  most  difficult  and  least 
measured  parameter  is  the  attentuation  coefficient  of 
radiation  as  it  penetrates  different  forest  canopies  (Nor- 
man and  Jarvis  1975,  Jarvis  et  al.  1976).  I  used  an  at- 
tenuation analysis  from  Kira  et  al.  (1969),  also  used  by 
}.  Rogers  (pers.  comm.)  in  his  general  watershed  hydro- 
logic  model  for  the  IBP  Coniferous  Forest  Biome.  Note 
that  his  analysis  attenuates  radiation  as  a  function  of 
LAI,  not  vertical  height,  eliminating  need  for  data  about 
canopy  geometry.  However,  implicit  in  this  treatment  is 
the  assumption  of  a  relatively  continuous  canopy.  A 
semiopen  stand  of  large  trees  may  have  deep  canopies 
but  an  overall  low  stand  LAI.  Also,  no  attempt  is  made  to 
adjust  radiation  effects  for  diurnal  sun-angle  changes  as 
is  done  in  more  rigorous  canopy  radiation  models  (Nor- 
man and  Jarvis  1975). 

G(12)  =  k,,  radiation  correction  fraction  (dim.) 

G(12)  =  1.0  -  (G(2)  -  G(ll))/G(2) 
IF(G(11)  .GT.  G(2))  G(12)  =  1.0 

where 

1.0     =  unity  multiplier  for  k,  reduction, 

G(2)    =  leaf  area  index, 

G(ll)  =  light-saturated  leaf  area  index, 

.GT.    =  a  FORTRAN  statement,  "greater  than." 

If  the  light-saturated  leaf  area  index  calculated  in 
G(ll)  is  larger  than  the  actual  LAI  of  the  stand,  then  no 
light-induced  k,  reduction  occurs,  and  G(12)  computes  a 
1.0  multiplier  for  k,.  If  light  is  limited  to  a  portion  of  the 
canopy,  G(12)  calculates  what  fraction  of  the  stand  LAI 
is  light-saturated  and  uses  that  fraction  as  a  multiplier 
for  the  k,  final  calculation.  Consequently,  if  70%  of  the 
canopy  is  receiving  0.1  ly  min-1  or  greater,  radiation  k, 
is  multiplied  by  0.7  to  reduce  the  canopy  average 
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because  of  the  canopy  that  has  closed  stomata.  Note 
that  this  assumes  a  k,  =  0  for  leaves  under  the  radiation 
threshold,  but  more  precise  estimates  of  k,  at  low-light 
levels  would  become  complicated  and  would  exceed  our 
knowledge  of  canopy  light  attenuation  theory  and  short- 
term  stomatal  fluctuations  under  low-light  conditions. 

G(13)  =  k,,  predawn  leaf  water  potential  correction 
(cm  sec  - ') 

G(13)  =  AMAX1  (B(12)  -  (B(12)/(B(14)  -  B(13)))-(G(10) 
-  B(13)),0.005) 
IF(G(10)  .LT.  B(13))  G(13)  =  B(12) 

where 

B(12)  =  maximum  canopy  average  k,, 
B(13)  =  spring  minimum  predawn  leaf  water  poten- 
tial, 

B(14)  =  stomatal  closure  threshold, 
G(10)  =  predawn  leaf  water  potential,  B^,, 
0.005  =  k,  at  complete  stomatal  closure, 
LT     =  a  FORTRAN  statement  "less  than." 
This  function  reduces  maximum  k,  as  soil  water  is 
depleted  (fig.  5).  Maximum  daily  k,  has  been  found  to  be 
a  direct  function  of  Bi/',,  G(10).  This  function  was  first 
demonstrated  on  Pseudotsuga  menziesii  in  Oregon  by 
Running  (1976)  and  duplicated  on  Pinus  contorta  in  Colo- 
rado by  Running  (1980),  both  on  sapling-size  trees. 
Although  I  hypothesize  that  the  relationship  between 
Bi/',  and  maximum  k,  is  fundamental,  different  species 
and  different  tree  sizes  may  require  a  modified  response 
curve.  Consequently,  B(12),  B(13),  and  B(14)  were  made 
externally  changeable.  The  ratio  B(12)/(B(14)-B(13))  pro- 
vides the  slope  of  the  k,  reduction  for  any  range  of  k,  and 
B(14)  entered. 

The  0.005  cm  s-1  minimum  k,  has  been  generally 
found  for  a  variety  of  conifers  (Hinckley  et  al.  1978).  It 
should  be  remembered  that  the  k,  value  calculated  here 
represents  a  canopy  average  for  leaves  of  different  ages 
through  all  crown  heights. 

Y  =  10.  x  (i-e-46(x)) 
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Q  I         ■         I         I         I  I  I  |  |  |  |  |  |  |  |  |  |  |  |  |  | 

0.0  0.5  1.0  1.5  2.0 

Incoming  shortwave  radiation  (ly  min"1  ) 

Figure  4.— Describes  the  leaf  area  index,  LAI,  that  would  receive  in- 
coming shortwave  radiation  in  excess  of  0.1  ly  min-1  (the 
stomatal  opening  threshold)  at  any  level  of  incoming  radiation. 


Y  =  0.15  -  0.136(X  ■  0.4) 


BY .  (MFa) 

Figure  5.— The  reduction  in  maximum  leaf  conductance,  k,, 
resulting  from  increasing  plant  water  stress  (predawn  leaf  water 
potential,  B^,). 

G(14)  =  k,,  humidity  correction  (cm  s-1) 

G(14)  =  G(13)  -  (G(13)- 0.05 -(G(l)-lO6- 4.0)) 
IF(G(14)  .LE.  0.0)  G(14)  =  0.005 

where 

G(13)  =  k,,  leaf  water  potential  correction, 

0.05    =  k,,  humidity  deficit  slope  coefficient  (fig.  6), 

G(l)    =  absolute  humidity  deficit, 

4.0     =  initial  humidity  reduction  threshold 

.LE.    =  a  FORTRAN  statement,  "less  than  or  equal 
to," 

0.005  =  minimum  kr 

There  is  much  evidence  that  k,  is  reduced  by  low  at- 
mospheric humidity,  regardless  of  \py  As  absolute 
humidity  deficit  calculated  in  G(l)  increases,  the  max- 
imum k,  from  G(13)  is  linearly  reduced  as  a  fraction  of 
the  original  (G(13))  leaf  conductance.  Theoretical 
analysis  of  this  phenomenon  is  reported  in  Lange  et  al. 
(1971),  Sheriff  (1977),  Edwards  and  Meidner  (1978),  and 
Losch  and  Schenk  (1978).  The  function  in  G(14)  is  based 
on  field  data  on  western  conifers  reported  by  Running 
(1976,  1980),  Kaufmann  (1976),  and  Tan  et  al.  (1977). 

From  reviewing  the  above  papers,  Running  (1980) 
found  that  this  function  reasonably  represents  field  data 
collected  on  Pseudotsuga  menziesii,  Pinus  ponderosa, 
Pinus  contorta,  Picea  sitchensis,  Picea  engelmannii,  and 
Tsuga  heterophyUa.  It  appears  that  when  humidity 
reduction  in  k,  is  analyzed  as  a  fraction  of  daily  max- 
imum kp  a  fairly  generalized  equation,  such  as  in  G(14), 
can  be  developed. 

Note  that  with  this  equation  structure  the  G(14)  equa- 
tion is  independent  of  maximum  k,,  which  is  input  from 
B(12)  through  G(13).  Also,  at  an  absolute  humidity  deficit 
above  25  ng  cm  3,  the  second  equation  sets  k,  = 
0.005  cm  s_1,  a  general  minimum  leaf  conductance  sub- 
stantiated by  many  studies.  Trees  in  very  arid  climates 
may  be  able  to  retain  stomatal  opening  at  these  humidity 
deficits. 

G(15)  =  k,,  air  temperature  correction  (cm  s_1) 
G(15)  =  S1(Z(3),  Z(7),G(14)) 

where 

Si  =  a  special  function  called  from  elsewhere  in 
program  that  uses  Z(3),  Z(7),  and  G(14)  values 
to  calculate  G(15), 
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Z(3)    =  air  temperature, 

Z(7)    =  hour  of  day, 

G(14)  =  k,,  humidity  reductions. 

This  equation  calls  special  function  Si  where  the  ef- 
fect of  air  temperature  of  k,  is  calculated  and  returned 
as  G(15). 

G(16)  =  kj,  final  canopy  average  (cm  s_1) 
G(16)  =  G(15)  •  G(12) 

where 
G(15)  =  k,, 

G(12)  =  radiation  limiting  fraction. 

This  computes  the  final  canopy  average  k,  by  the 
radiation  reduction  fraction  from  G(12). 

To  summarize  k,  control  in  the  model:  G(13)  computes 
leaf  (and  soil)  water  potential  effects.  G(14)  adds  humid- 
ity effects  to  kr  G(15)  adds  air  temperature  effects. 
Finally,  G(16)  adds  radiation  control  from  G(12).  The 
overall  hierarchy  of  k,  control  from  strongest  to  weakest 
is  (1)  radiation,  (2)  leaf  water  potential  or  B^,,  (3) 
humidity  or  ABSHD,  and  (4)  air  temperature. 

G(17)  =  transpiration  flux  density  (g  cm  2  s_1) 
G(17)  =  G(16)  •  G(l) 

where 

G(16)  =  kj,  canopy  average, 

G(l)    =  absolute  humidity  deficit  (ABSHD). 

Transpiration  flux  density  is  the  product  of  k,  and  AB- 
SHD. This  equation  is  equivalent  to  the  diffusion  equation 
used  by  Tan  et  al.  (1978). 

E  =  (PCp/LV)((e  -  ea)/(rs  +  rj 

where 

E    =  transpiration  rate  (g  cm"2  s_1), 

P     =  air  density  (g  cm-3), 

Cp  =  specific  heat  of  air  Q  g_1  °C  ^, 

L    =  latent  heat  of  vaporization  Q  g_1), 

V    =  psychrometric  constant  (mb  °C 

es    =  saturation  vapor  pressure  (mb), 

ea    =  ambient  vapor  pressure  (mb), 

rs    =  stomatal  resistance  (s  cm 

ra    =  aerodynamic  resistance  (s  cm"1). 


ABDHO 

0  20r 


ABSHD  (g  cm  3  x  106) 

Figure  6.— The  reduction  in  maximum  leaf  conductance,  k 
resulting  from  increasing  absolute  humidity  deficit. 


In  the  equation  in  G(17),  P,  Cp,  L,  V,  es,  and  ea  are  all 
subsumed  into  G(l),  the  absolute  humidity  deficit.  Both 
Tan  et  al.  (1978)  and  the  research  reported  here  ignored 
ra  for  the  reason  that  ra  is  typically  20-100  times  smaller 
than  rs  in  coniferous  forests  (Smith  1980),  and  would  re- 
quire input  of  vertical  profiles  of  windspeed.  Also,  k, 
=  1/r  is  used. 

s 

Midday  humidity  in  western  forests  is  so  low  that  the 
humidity  deficit  component  of  the  Penman  -  Monteith 
equation  far  exceeds  the  radiation  component  of  evapo- 
rative demand  (Jarvis  et  al.  1976,  Luxmoore  et  al.  1981). 
Consequently,  the  diffusion  equation  used  in  G(17)  to 
compute  stand  transpiration  and  used  by  Tan  et  al. 
(1978)  in  their  transpiration  model  varies  from  the  more 
theoretically  satisfying  Penman  -  Monteith  equation  by 
less  than  2%.  The  diffusion  equation  is  used  because  of 
its  simplicity  and  reduced  demand  for  meteorological  in- 
put variables. 

G(18)  =  hourly  stand  transpiration  (g  h-1) 
G(18)  =  G(17)-B(9)-3600 

where 

G(17)  •    transpiration  flux  density, 
B(9)    =  stand  leaf  area, 
3600  =  seconds  per  hour. 

G(18)  equals  F(8,9),  water  transfer  from  deaf  storage, 
X(8),  to  the  atmospheric  transpiration  sink,  X(9). 


Soil  Water  Input 

G(50)  =  effective  precipitation  (cm3  h_1) 


IF(Z(7)  .GT.  4.)  Z(2)  =  0.0 


G(50)  = 

=  AMAXl((Z(2)/4-B(6)  ■  (B(9)/B(10)))  ■  B(10),0.0) 

where 

Z(7)  = 

hour  of  day, 

Z(2)  = 

precipitation, 

B(6)  = 

canopy  interception  coefficient, 

B(9)  = 

stand  leaf  area, 

B(10)  = 

ground  surface  area. 

This  function  translates  the  one-dimensional  precipi- 
tation data  into  volumetric  units  and  subtracts  canopy 
interception  based  on  leaf  area  index,  B(9)/B(10).  The 
result  is  equal  to  F(3,3),  the  precipitation  input  into  X(3). 

Precipitation  Z(2)  is  divided  by  4  to  estimate  input 
from  the  once-daily  precipitation  record  and  is  entered 
between  0000  and  0400  hours.  This  may  underestimate 
instantaneous  rainfall  intensity  and  overestimate  in- 
terception losses.  With  modern  data-loggers,  hourly 
precipitation  should  be  recorded  to  avoid  these  artificial 
problems.  Canopy  interception  of  precipitation  is 
treated  proportionally  to  LAI  using  an  interception  coef- 
ficient which,  in  effect,  is  leaf  surface  water  storage. 
Nearly  two-thirds  of  the  annual  precipitation  in  a  Rocky 
Mountain  lodgepole  pine  forest  falls  as  snow,  outside  of 
the  growing  season  (Alexander  and  Watkins  1977).  Sum- 
mer rainfall  typically  is  in  the  form  of  short  thunder- 
showers  that  do  not  recharge  the  rooting  zone  soil  and 
are  a  modest  component  of  the  stand  hydrologic 
balance.  Consequently,  while  more  elaborate  canopy  in- 
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terception  and  evaporation  models  are  available  (Rutter 
et  al.  1971,  Murphy  and  Knoerr  1975,  Stewart  1977, 
Ford  and  Deans  1978),  a  more  accurate  treatment  of  in- 
terception processes  does  not  seem  warranted  here. 
This  algorithm  assumes  that  the  canopy  is  dry  at  each 
hourly  iteration,  which  overestimates  interception  loss 
during  an  extended  rainfall  event.  Also,  no  attempt  is 
made  to  retard  transpiration  while  the  canopy  is  wet. 
However,  this  problem  is  reduced  implicitly  because, 
after  a  thundershower,  evaporative  demand  is  greatly 
reduced  until  the  intercepted  water  has  evaporated. 

Application  of  this  model  to  sites  having  more  signifi- 
cant summer  precipitation  should  entail  collecting  hour- 
ly precipitation  data  and  possibly  calculating  a  canopy 
energy  balance  to  predict  evaporation  of  intercepted 
precipitation. 

The  AMAXl  function  with  a  0.0  limit  is  used  to  protect 
the  equation  from  calculating  a  negative  precipitation. 

G(51)  =  snowmelt  (cm3  h_1) 

G(51)  =  AMIN1(Z(3)  •  B(l)  •  B(10),  X(l)) 
IF(Z(3)  .LE.  0.0)  G(51)  =  0.0 

where 

Z(3)    =  air  temperature, 

B(l)    =  snowpack  melt  coefficient, 

B(10)  =  ground  surface  area, 

X(l)    =  snowpack  water  content. 

This  equation  melts  the  snowpack  in  X(l)  as  a  function 
of  air  temperature,  Z(3),  and  a  melt  coefficient,  B(l), 
then  multiplies  it  by  ground  area,  B(10),  to  give 
volumetric  water  input  to  X(3).  The  degree-day  snow- 
melt  equation  was  basically  derived  from  data  in  the 
U.S.  Army  Corps  of  Engineers  Snow  Hydrology  report 
(1956).  This  approach  is  only  reasonable  for  spring 
snowmelt  of  a  ripe  isothermal  pack.  More  rigorous 
models  of  snowmelt  do  complete  energy  balances  for  the 
snowpack  (Leaf  and  Brink  1973,  Colbeck  and  Ray  1979). 
This  increased  complexity  seems  unnecessary  for  a 
model  whose  focus  is  primarily  summer  water  stress 
development.  The  melt  coefficient,  B(l),  should  be  recali- 
brated for  use  in  other  regions.  G(51)  is  equal  to  F(l,3). 

G(52)  =  litter  vapor  conductance  (cm  s-1) 
G(52)  =  EXP(5.0  •  (G(3)  -  1.0)) 

where 

5.0,  1.0  =  curve-fitting  coefficients, 
G(3)      =  litter  water  fraction. 

This  equation  calculates  a  water  vapor  conductance 
term  for  evaporation  of  water  from  the  forest  floor  litter. 
My  approach  is  derived  from  classic  evaporation  and 
diffusion  principles  (Lee  1980),  and  is  similar  to  the  litter 
evaporation  function  of  CFBMG  (1979),  except  that  a 
variable  conductance  term,  G(52),  was  calculated  as  a 
function  of  litter  water  content.  It  assumes  a  conduct- 
ance of  1.0  cm  s_1  from  a  saturated  litter  layer,  with  an 
exponential  decrease  in  conductance  as  the  litter  sur- 
face dries,  increasing  boundary  layer  resistance  (fig.  7). 

G(53)  =  litter  evaporation  (cm3  h_1) 

G(53)  =  G(52)-G(l)-B(10)-3600 


o1      1      1  — 1  1  ■  1  ■  1  ■  ' 

0.0  0.2  0  4  06  0.8  1.0 

Forest  floor  litter  water  content  (%  of  capacity) 


Figure  7.— As  water  content  of  the  forest  floor  titter  is  depleted,  the 
water  vapor  conductance  is  reduced  by  this  function. 

IF(G(3)  .LE.  0.0)  G(53)  =  0.0  , 
IF(X(1)  .GT.  0.0)  G(53)  =  0.0 

where 

G(52)  =  litter  vapor  conductance, 

G(l)    =  absolute  humidity  deficit, 

B(10)  =  ground  surface  area, 

3600  =  seconds  per  hour, 

G(3)    =  litter  water  fraction, 

X(l)    =  snowpack  water  content. 

Litter  evaporation  is  driven  by  atmospheric  demand, 
absolute  humidity  deficit,  G(l),  and  is  limited  by  G(52)  as 
the  surface  dries.  Volumetric  litter  evaporation  per  hour 
is  determined  by  multiplying  by  B(10)  and  3,600  s/hour. 
This  approach  is  simple  and  can  be  calibrated  easily. 
Model  adjustments  would  best  be  made  in  the  G(52) 
response  curve  of  figure  7.  Because  this  equation  only 
quantifies  mass  transfer  considerations  of  evaporation 
theory,  water  loss  from  an  open  surface  subjected  to  a 
significant  diurnal  radiation  load  may  be  underesti- 
mated. G(53)  is  subtracted  from  the  water  input  to  X(3) 
by  F(3,3). 

G(54)  =  surface  runoff  (cm3  h_1) 

G(54)  =  AMAXl((G(50)  +  G(51)-(B(2)-B(10))),  0.0) 
where 
G(50)  =  precipitation, 
G(51)  =  snowmelt, 

B(2)    =  maximum  soil  surface  infiltration  rate, 
B(10)  =  ground  surface  area. 

If  precipitation,  G(50),  plus  snowmelt,  G(51),  exceed 
the  soil  surface  infiltration  capacity,  B(2)-B(10),  runoff 
occurs  to  X(2).  This  water  transfer  is  done  in  F(2,2). 
However,  overland  flow  on  forested  land  is  rare  (Ander- 
son et  al.  1976).  In  addition,  unless  precipitation  data  is 
input  hourly  instead  of  the  present  daily  data,  there  is 
little  potential  for  overland  flow  in  the  model. 

G(55)  =  litter  compartment  water  transfer  (cm3  h_1) 

G(55)  =  AMAX1((G(50)  +  G(51)-G(53) 
-G(54))-(B(3)-X(3)),  0.0) 
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where 

G(50)  = 

precipitation, 

G(51)  = 

snowmelt, 

G(53)  = 

litter  surface  evaporation, 

G(54)  = 

surface  runoff, 

B(3)  = 

water  storage  capacity  of  forest  floor  litter, 

X(3)  = 

current  water  content  of  forest  floor  litter. 

This  function  does  a  mass  balance  of  the  forest  floor 
litter  storage,  X(3),  to  determine  the  volume  of  water 
flow  from  X(3)  into  X(4).  Inputs  of  precipitation,  G(50), 
and  snowmelt,  G(51),  are  added;  outputs  of  litter  evapo- 
ration, G(53),  and  surface  runoff,  G(54),  are  subtracted. 
If  X(3)  is  not  at  capacity,  the  deficit,  B(3)-X(3),  is  first 
refilled,  with  the  excess  input  water  then  cascading  into 
X(4).  The  result  equals  F(3,4). 

G(56)  =  root  zone  Soil  A  water  transfer  (cm3  h_1) 
G(56)  =  AMAX1(G(55)-(B(4)-X(4)),0.0) 

where 

G(55)  =  litter  compartment  water  transfer, 

B(4)    =  maximum  water  storage  in  root  zone  Soil  A, 

X(4)    =  current  water  storage  in  root  zone  Soil  A. 

As  in  G(55),  a  mass  balance  is  done  for  X(4).  Input 
from  X(3)  is  added  as  G(55),  water  used  to  refill  any 
deficit,  B(4)-X(4),  is  subtracted,  and  the  remainder 
cascaded  to  X(5).  G(56)  equals  F(4,5). 

G(57)  =  root  zone  Soil  B  water  transfer  (cm3  h_1) 
G(57)  =  AMAX1(G(56)-(B(5)-X(5)),0.0) 

where 

G(56)  =  root  zone  Soil  A  water  transfer, 

B(5)    =  maximum  water  storage  in  root  zone  Soil  B, 

X(5)    =  current  water  storage  in  root  zone  Soil  B. 

Again,  input  from  X(4)  is  added  until  any  deficit  in  X(5) 
is  filled,  B(5)-X(5).  The  excess  water  cascades  into  the 
subsurface  outflow  compartment,  X(6),  with  an  un- 
limited capacity.  G(57)  equals  F(5,6). 

Inherent  in  the  mass  balances  of  G(55),  G(56),  and  G(57) 
are  two  assumptions.  First,  they  do  not  allow  soil  super- 
saturation.  All  water  in  excess  of  field  capacity  is  routed 
to  the  lower  compartment.  A  second  related  assumption 
is  that  they  effectively  treat  saturated  hydraulic  conduc- 
tivity of  the  soil  as  infinite.  At  each  hourly  iteration,  ex- 
cess water  is  immediately  routed  to  the  lower  compart- 
ment with  no  lag  time.  This  was  done  to  reduce  the  model 
requirements  for  soil  physical  data,  and  because  these 
rather  short-term  phenomena  were  of  small  importance 
to  seasonal  water  stress  development.  In  some  areas  and 
for  certain  applications,  this  part  of  the  model  may  re- 
quire more  sophisticated  treatment. 

G(58)  =  canopy  evaporation  (cm3  h1) 

G(58)  =  Z(2)/4.-B(10)-G(50) 

where 

Z(2)/4.  =  precipitation,  entered  from  0000  to  0400 
hours, 

B(10)    =  ground  surface  area, 
G(50)    =  effective  precipitation. 


Water  Withdrawal 

G(60)  =  leaf  water  storage,  current  available  (cm3) 
G(60)  =  B(8)-(l.-G(10)/B(14)  +  B(13)/B(14)) 

where 

B(8)    =  maximum  available  leaf  storage, 

G(10)  =  predawn  leaf  water  potential,  B^,, 

B(14)  =  leaf    osmotic    potential,    stomatal  closure 

threshold, 
B(13)  =  spring  minimum  B\pv 

Recent  research  on  the  role  of  internal  water  storage 
in  plant  water  relations  (Hellkvist  et  al.  1974,  Roberts 
1977,  Richter  1978,  Waring  and  Running  1978,  Waring 
et  al.  1979,  Running  1980c)  has  demonstrated  that  leaf 
water  content  and  leaf  water  potential  are  directly 
related.  This  equation  reduces  available  leaf  water 
storage  as  predawn  ^  increases,  G(10),  and  approaches 
the  leaf  osmotic  potential,  B(14),  where  stomata  close. 

G(61)  =  leaf  storage  exchange  rate  (cm3  h-1) 
G(61)  =  X(8)/B(16) 

where 

X(8)    =  leaf  water  storage, 

B(16)  =  leaf  water  depletion  coefficient. 

Although  available  leaf  water  storage  has  been 
studied  in  some  detail,  the  short-term  exchange  rates  of 
water  from  leaf  tissue  has  not.  The  flow  of  water 
through  leaf  tissue  greatly  exceeds  the  diurnal  net 
change  in  leaf  water  content.  This  function,  taken  from 
Running  (1980c)  on  Pinus  contorta,  limits  the  hourly  ex- 
change rate  of  water  from  leaf  storage,  X(8),  and  is 
necessary  for  both  biological  and  mathematical  stability 
of  leaf  water  exchange. 

G(62)  =  stem  water  transport  (cm3  h"1) 

G(62)  =  0.0 
IF(G(61)  .GE.  G(18))  G(62)  =  AMAX1(G(60) 
-X(8),G(18)/B(16)) 
IF(G(61)  .LT.  G(18))  G(62)  =  G(18)-G(61) 

where 

G(61)  =  leaf  storage  exchange  rate, 

G(18)  =  hourly  stand  transpiration, 

G(60)  =  current  leaf  water  storage, 

X(8)    =  leaf  water  storage, 

B(16)  =  leaf  water  depletion  coefficient. 

This  function  computes  the  water  flow  from  X(7)  to 
X(8)  by  determining  how  much  of  the  transpiration  de- 
mand can  be  satisfied  by  X(8),  because  in  order  to  retain 
mass  balance  in  the  model,  all  transpiration  demands 
must  be  met.  When  demand  generates  water  stress  by 
water  depletion  in  various  parts  of  the  model,  transpira- 
tion is  reduced  in  future  iterations  by  feedback  controls 
on  k,.  If  exchangeable  water  in  X(8),  G(61),  exceeds  the 
transpiration  demand  G(18),  then  G(62)  can  either  refill 
previous  X(8)  deficits  G(60)-X(8)  or  a  fraction  of  the  de- 
mand G(18)/B(16).  In  the  more  normal  case,  where  de- 
mand G(18)  exceeds  supply  G(61),  then  G(62)  requests 
flow  from  X(7)  to  cover  all  demand  above  exchangeable 
supply  in  X(8),  G(18)-G(61).  G(62)  equals  F(7,8). 
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G(66)  =  stem  storage  exchange  rate  (cm3  h  ') 
G(66)  =  X(7)/B(17) 

where 

X(7)    =  stem  water  storage, 

B(17)  =  stem  water  exchange  coefficient. 

From  research  by  Waring  and  Running  (1978),  War- 
ing et  al.  (1979),  and  Running  (1980c),  there  is  data  on 
three  different  coniferous  species  showing  sapwood 
water  depletion  rates.  This  function  limits  the  amount  of 
water  that  can  be  withdrawn  from  X(7)  on  any  one  itera- 
tion by  a  rate  constant  B(17). 

G(67)  =  root  zone  Soil  A  water  uptake  (cm3  h  ') 

G(67)  =  AMIN1(G(62)-G(66),  (B(7)-X(7)/B(17),  X(4)/B(17)) 
where 

G(62)  =  stem  water  transport, 

G(66)  =  stem  storage  exchange  rate, 

B(7)  =  maximum  available  stem  water  storage, 

X(7)  =  current  stem  water  storage, 

B(17)  =  stem  and  soil  water  exchange  coefficient, 

X(4)  =  root  zone  Soil  A  soil  water. 

Flow  from  X(4)  into  X(7)  is  determined  by  the  demand, 
G(62),  minus  any  that  can  be  satisfied  by  X(7),  G(66).  If 
there  is  no  demand  (i.e.,  at  night),  but  X(7)  is  not  full,  flow 
will  come  into  X(7)  to  fill  this  deficit,  B(7)-X(7),  subject  to 
the  B(17)  exchange  coefficient.  The  final  control  of  F(4,7) 
uptake  is  that  water  withdrawal  cannot  exceed  the 
exchangeable  supply  in  rooting  zone  Soil  A,  X(4)/B(17). 
G(67)  equals  F(4,7). 

G(68)  =  root  zone  Soil  B  water  uptake  (cm3  h  ') 


G(68)  = 

AMIN1(G(62)-G(66)-G(67),  (B(7)-X(7))/B(17)) 

where 

G(62)  = 

stem  water  transport, 

G{66)  = 

stem  storage  exchange  rate, 

G(67)  = 

root  zone  Soil  A  water  uptake, 

B(7)  = 

maximum  available  stem  water  storage, 

X(7)  = 

current  stem  water  storage, 

B(17)  = 

stem  water  exchange  coefficient. 

Similar  to  G(67),  G(68)  is  the  water  requirement 
needed  to  fulfill  transpiration  demand  minus  supply  al- 
ready contributed  by  stem  water,  G(66),  and  by  X{4), 
G(67).  This  function  chooses  the  minimum  of  two  calcu- 
lated flow  potentials — the  remaining  demand  or  the 
exchangeable  supply  in  X(7).  If  X(7)  is  not  full,  the  deficit 
B(7)-X(7)  can  also  be  filled  subject  to  the  exchange  coef- 
ficient limit  B(17). 

Water  uptake  is  controlled  by  a  feedback  loop  where 
soil  moisture  depletion,  X(4)  and  X(5),  increase  predawn 
plant  moisture  stress,  G(10),  which  reduces  k,,  G(13),  and 
transpiration,  G(18),  thus  reducing  water  uptake  de- 
mand, G(67)  and  G(68). 

If  transpiration  demand  exceeds  supply  in  X(4),  X(7), 
and  X(5),  which  could  only  happen  after  months  of  pro- 
gressive tree  water  stress  and  depletion  of  soil  moisture 
to  zero,  X(5)  would  become  negative.  This  means  the  tree 
or  stand  died  of  water  stress! 


G(69)  =  total  root  system  water  uptake  (cm3  h" ') 
G(69)  =  G(67)  +  G(68) 

where 

G(67)  =  root  zone  Soil  A  water  uptake, 
G(68)  =  root  zone  Soil  B  water  uptake, 
Total  root  system  uptake  is  combined  for  use  in 
calculating  diurnal  leaf  water  potential  in  G(73). 

Leaf  Water  Potential,  Flow  Resistances,  Transpiration 
Ratios 

G(70)  =  root  resistance,  soil  temperature  influences 
(MPa  /xg  1  cm  2  s_1) 

IF(Z(5)  .LE.  0.0)  G(70)  -  10.0 
IF(Z(5)  .GT.  0.0)  G(70)  =  1.0/Z(5) 

where 

Z(5)  =  soil  temperature, 

10.0  =  a  high  soil  resistance  limit  when  soil 
temperature  is  less  than  0°  C;  also  protects 
the  fraction  from  division  by  0.0. 

Recent  studies  have  found  root  resistance  to  water 
uptake  to  increase  rapidly  as  soil  temperature  ap- 
proaches freezing  (Babalola  et  al.  1968,  Havranek  1972, 
Kaufmann  1975,  Running  and  Reid  1980).  These  effects 
have  been  found  in  different  conifer  species  and  are 
partially  physiological  in  nature.  Increasing  water 
viscosity  explains  only  a  part  of  the  observed  increase  in 
flow  resistance.  The  function  in  figure  8  approximates 
the  data  of  Running  and  Reid  (1980)  on  Pinus  contorta 
and  of  Kaufmann  (1975)  on  Picea  engelmannii. 

G(71)  =  flow  resistance,  soil  water  influences  (MPa  /xg_1 
cm-2  s_1) 

G(71)  =  -  0.28  +  2.0  •  G(10) 

where 

-0.28,  2.0  =  empirical  coefficients  to  fit  data  from 
Running  (1980b), 

G(10)        =  predawn  leaf  water  potential. 

The  importance  of  flow  resistance  in  the  overall  devel- 
opment of  water  potential  gradients  in  coniferous  trees 
has  been  quantified  using  a  root  excision  technique 
(Roberts  1977,  Running  1980b).  Root  resistance  appears  to 
increase  as  a  direct  function  of  decreasing  soil  water 
potential.  This  function  increases  flow  resistance  as 
predawn  leaf  water  potential  increases  (negative  \p  values 
are  treated  as  positive  numbers  in  the  model,  fig.  9).  These 
data  were  from  trees  3-5  m  tall.  Both  the  relaive  and  ab- 
solute influence  of  root  resistance  probably  will  be  dif- 
ferent on  larger  trees. 

G(72)  =  total  root  resistance  (MPa  /xg~'  cm~2  s_1) 
G(72)  =  G(70)  +  G(71) 

where 

G(70)  =  root  resistance,  soil  temperature  influences, 
G(71)  =  root  resistance,  soil  water  influences. 

G(73)  =  leaf  water  potential,  ^  (MPa) 
G(73)  =  AMINl(G(10)  +  G(72)-(G(69)-278.0/B(9)),  B(14)) 
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where 

G(10)  = 

predawn  leaf  water  potential, 

G(72)  = 

total  root  resistance, 

G(69)  = 

total  root  water  uptake, 

278.0  = 

combines  the  constants  10Vg  g_1  and 

3,600  s  h-\ 

B(9)  = 

stand  leaf  area, 

B(14)  = 

leaf   osmotic    potential,    stomatal  closure 

threshold. 

From  plant  water  relations  theory: 

T  =  MsoU  -  A  eaf)/Rtotal 

where 

T 

=  transpiration  rate, 

^soil  ^leaf 

=  the  soil  to  leaf  water  potential  gradient, 

B-total 

=  the  total  resistance  to  water  flow  between 

the  soil  and  leaf. 

This  equation  can  be  rearranged  to  predict  i/*,  ,  as 
follows:  \j/]eaf  -  ^soi]  +  Rtotal'T  [\p  treated  as  positive  value) 
and  has  been  used  in  a  number  of  studies  (Elfving  et  al. 
1972,  Kaufmann  and  Hall  1974,  Hinckley  et  al.  1978). 
However,  as  pointed  out  by  Waring  and  Running  (1976, 
1978)  and  Running  (1980c),  system  capacitance  confuses 
this  equation.  I  have  modified  the  equation  by 
substituting  G(10)  for  ^soU  and  by  using  G(69)  instead  of 
transpiration  G(18)  for  water  flux.  This  allows  water  to  be 
withdrawn  from  internal  reservoirs  before  \p{  begins  to  in- 
crease and  produces  the  diurnal  ^  hysteresis  typically 
found  in  trees  (Jarvis  1975,  Running  1980a).  Using  root 
water  uptake,  G(69),  to  generate  \f/}  also  agrees  with  studies 
showing  the  major  component  of  overall  flow  resistance, 
and  \J/  gradient  development  to  be  at  the  soil-root  interface 
(Robert  1977,  Running  1980b). 

Development  of  \J/{  is  limited  in  this  equation  to  B(14), 
the  stomatal  closure  threshold,  which  typically  varies  bet- 
ween 1.5  and  2.3  MPa  for  temperate  conifers  (Hinckley  et 
al.  1978).  Under  severe  stress,  trees  can  exceed  this  value, 
although  the  model  cannot  presently  handle  that 
circumstance. 


.2 1  i  i  i  i  I  i  i  i  i  I  i  i  i  i  I  i  i  i  i  I  i  i  i  i_J 

0  1.0  2.0  3.0  4.0  5.0 

BY,  (MPa) 

Figure  9.— As  tree  water  stress,  B</<,,  increases,  the  total  water  flow 
resistance  in  the  tree  increases. 

where 

G(16)  =  k,,  canopy  average, 
B(12)  =  maximum  average  kr 

The  T/PT  ratio  provides  a  convenient  quick  measure  of 
the  current  degree  of  stomatal  closure  controlling  water 
loss  (Reed  and  Waring  1974).  This  equation  simplifies  the 
T/PT  ratio  by  removing  the  humidity  deficit  term  from  the 
numerator  and  denominator. 

G(77)  =  areal  transpiration  (mm  h_1) 

G(77)  =  G(18)/B(10)-10 

where 

G(18)  =  stand  transpiration, 
B(10)  =  ground  surface  area, 

10      =  conversion  factor  from  cm3  h1  to  mm  h_1. 
This  calculates  transpiration  in  simple  water  depth 
equivalence,  as  is  routinely  used  in  hydrology  literature. 


transpiration/potential  transpiration  ratio  (dim.) 
G(76)  =  G(16)/B(12) 


1/X 


SPECIAL  FUNCTIONS 

Special  Function  Si — Frost  k.  Function 

Function  SlfTAIR,  XHOUR,  COND) 
IF(XHOUR  -  5.0)  91,  90,  91 

90  SRT  =  TAIR 

91  IF(SRT  .  GT  .  0.0)  GO  TO  92 

Si  =  AMAXl(COND  +  0.02  •  SRT,  0.005) 
IF(SRT  .  LT  .  -  7.0)  Si  =  0.005 
GO  TO  93 

92  Si  =  COND  +  COND  •  0.003  •  (TAIR  -  10.0) 

93  RETURN 
END 


G(76) 
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Figure  8.— Decreasing  soil  temperature  produces  a  nonlinear  in- 
crease in  root  water  flow  resistance. 


The  effect  of  air  temperature  on  k,  is  calculated  in 
Function  Si  and  called  by  G(15).  The  temperature  effect 
has  two  sections.  When  air  temperature  is  above  10°  C, 
k,  is  increased  modestly  (0.003  cm  s_1  °C_1),  if  below 
10°  C,  k,  is  decreased  by  the  same  coefficient  (Hinckley 
et  al.  1978).  There  is  evidence  that  subfreezing  night 
temperatures  retard  stomatal  opening  even  after  tem- 
perature recovery  (Neilson  and  Jarvis  1976,  Fahey  1979, 
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Kaufmann  1982);  therefore,  the  other  part  of  Function 
Si  determines  if  frost  occurred  at  dawn,  defined  as 
0500  hours  (SRT).  If  it  had  frozen,  a  much  stronger  k, 
reduction  is  used  (0.02  cm  s_1  °C-')  until  0500  the  next 
morning. 

Special  Function  S3 — Absolute  Humidity  Deficit 
Function  S3  (TA,  RH) 

ESD  =  6.1078  •  EXP((17.269  ■  TA)/(237.3  +  TA)) 

ES     =  RH/100  •  ESD 

VPD  =  AMAX1((ESD  -  ES),  0.0) 

S3      =  217.0  E  -6  •  VPD/(TA  +  273.16) 

RETURN 

END 

If  the  evaporating  surface  of  a  leaf  is  assumed  to  be  at 
saturation  and  leaf  temperature  to  equal  air  tempera- 
ture, the  evaporative  demand  in  mass  units  can  be  ob- 
tained by  first  calculating  saturation  vapor  pressure  at 
TAIR  using  the  equation  of  Murray  (1967).  Then  calcu- 
late ambient  vapor  pressure  (ES)  using  input  relative 
humidity  (RH).  Vapor  pressure  deficit  is  defined  as  (ESD 
-  ES).  Absolute  humidity  deficity  (ABSHD)  in  grams  per 
cubic  centimeter  is  then  calculated  from  VPD  and  TA 
and  returned  to  G(l)  as  S3.  This  function  can  be  modified 
easily  to  calculate  ABSHD  from  dewpoint  or  wetbulb 
depression  instead  of  relative  humidity. 

DAYTRANS 

The  primary  purpose  of  DAYTRANS  is  to  develop  a 
more  simplified,  less  data-intensive  model  that  retains  as 
much  of  the  prediction  accuracy  of  H20TRANS  as  possi- 
ble. A  major  simplification  was  going  to  a  daily  timestep 
in  DAYTRANS.  As  a  result,  computer  execution  cost  is 
also  greatly  reduced.  The  overall  structure  and  perspec- 
tive of  the  models  are  similar.  DAYTRANS  was  written 
after  completion  of  H20TRANS  by  removing  the  less  crit- 
ical sections  of  H20TRANS.  Because  much  of  the  docu- 
mentation of  DAYTRANS  duplicates  H20TRANS,  discus- 
sion here  will  concentrate  on  the  differences  between  the 
models.  Readers  can  refer  to  the  comprehensive 
H20TRANS  section  for  documentation  of  the  sections 
that  are  basically  identical.  See  figure  2  for  the 
DAYTRANS  compartment  flow  diagram. 


STATE  VARIABLES— X(i) 

A  list  of  state  variables  follows.  See  definitions  for  the 
comparable  state  variables  in  the  H20TRANS  section. 


Unit 

X(l)  = 

snowpack  water  content 

cm3 

X(2)  - 

root  zone  soil  water  content 

cm3 

X(3)  = 

tree  sapwood  water  content 

cm3 

X(4)  = 

transpiration  (unlimited  capacity) 

cm3 

X(5)  = 

surface  runoff  (unlimited  capacity) 

cm3 

X(6)  = 

subsurface  outflow  (unlimited  capacity) 

cm3 

X(7)  = 

evaporation  (unlimited  capacity) 

cm3 

FLOW  FUNCTIONS— F(i,j) 

A  list  of  flow  functions  follows.  See  F(i,j)  section  in 
H20TRANS  for  explanations. 
F(l,2)  =  snowmelt  input,  G(50) 
F(2,2)  =  precipitation  input,  G(l) 
F(2,3)  =  root  water  uptake,  G(53) 
F(2,5)  =  surface  runoff,  G(51) 
F(2,6)  =  subsurface  outflow,  G(52) 
F(3,4)  =  transpiration,  G(19) 
F(7,7)  =  evaporation,  canopy  and  litter,  G(9) 


DRIVING  VARIABLES— Z(i) 

DAYTRANS  uses  the  following  driving  variables. 


Unit 

Z(l)  = 

Julian  date 

day 

Z(2)  = 

precipitation 

cm 

Z{3)  = 

air  temperature 

(daylight  average) 

°C 

Z{4)  = 

absolute  humidity  deficit 

(daylight  average) 

g  cm_; 

Z(5)  = 

soil  temperature,  20  cm 

°C 

Z(6)  = 

incoming  shortwave  radiation 

(daylight  average) 

ly  min 

Z(7)  = 

daylength  (sunrise-sunset) 

s 

Z(8)  = 

night  minimum  air 

temperature 

°C 

A  primary  reason  for  using  DAYTRANS  occurs  when 
only  daily  meteorological  data  is  available.  It  became 
evident  executing  both  models  that  a  critical  factor  in 
the  success  of  the  DAYTRANS  simulation  was  the 
accuracy  with  which  the  climatic  variables  were  aver- 
aged. As  more  meteorological  stations  record  on  micro- 
processor data-loggers,  daily  averages  can  be  gener- 
ated by  electronic  integration  of  the  variable  inputs 
throughout  the  day.  Some  suggestions  on  averaging 
certain  daily  meteorological  data  for  use  in  DAY- 
TRANS are  presented  next.  Note  that  Z(3),  Z(4),  and 
Z(6)  must  be  averaged  for  the  daylight  period  only.  This 
is  based  on  the  assumption  that  transpiration  and 
water  stress  development  occur  only  during  the  day. 

Air  Temperature — Z(3) 

Frequently,  only  maximum  and  minimum  air  tempera- 
tures are  recorded  for  a  day.  If  the  daylight  course  of  air 
temperature  is  assumed  to  approximate  three  quadrants 
of  a  sine  function  (Parton  and  Logan  1981)  as  shown  in 
figure  10,  a  simple  approximation  of  daylight  average  air 
temperature  can  be  derived. 
The  arithmetic  mean  is: 

T  =  (T     +  T.J/2  =  sin  0. 

v    max  miiv 

But  integrating  the  sine  function  yields 

PI 

1  (      sin  x  dx  =  0.212 

PI -(-PI/2)  _J 

where 
PI  =  3.1416. 
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So 

T     =  0.212(T     -  T)  +  T 

ave  v    max  ' 

where 

Tave  =  average  daylight  air  temperature, 
T     =  arithmetic  mean  of  T     and  T  . 

max  mm 

This  equation  was  used  to  generate  a  DAYTRANS 
data  file  that  was  compared  against  a  data  file 
developed  from  averaging  recorded  hourly  air  tempera- 
ture values.  The  regression  of  hourly  averaged  air 
temperature  daylight  averages  against  the  sine  function 
weighted  daylight  averages  was: 

Tha  =  -1.14  +  1.12  T  f   R2  =  0.93  n  =  120  days 
where 

Tha   =  air  temperature,  hourly  averaged, 
T  f    =  air  temperature,  weighted  sine  function. 
The   same   regression   against   arithmetic  maximum- 
minimum  averages  was 

T    =  0.74  +  1.20  T       R2  =  0.88  n  =  120  days 

ha  mm  J 

where 

Tmm  =  maximum-minimum  averaged  temperatures. 

The  sine  function-weighted  average  more  closely  ap- 
proaches the  optimum  Tave  =  Tha  prediction,  as  shown 
by  a  comparison  of  Y-intercepts,  slopes  and  correlation 
coefficients. 

It  should  be  emphasized  that  these  data  manipula- 
tions are  done  in  developing  the  DAYTRANS  input  data 
file.  This  is  not  done  internally  by  the  DAYTRANS 
program,  although  it  could  be  if  a  permanent  data  base 
required  it. 

Absolute  Humidity  Deficit— Z(4) 

Absolute  humidity  deficit  is  calculated  for  Z(4)  by 
reading  in  TAIR,  air  temperature,  and  XRHUM  (relative 
humidity).  Four  equations  are  then  used  to  calculate 
saturation  vapor  pressure,  absolute  humidity,  and  ab- 
solute humidity  deficit.  These  equations  are  covered  in 
Special  Function  S3  in  H20TRANS.  Some  considerations 
for  estimating  daily  average  relative  humidity  follow. 

The  diurnal  course  of  relative  humidity  often  is  a  mir- 
ror image  of  the  air  temperature  trace  because  of  the 
influence  of  air  temperature  on  saturation  vapor  pres- 
sure. Consequently,  a  sine  function-weighted  correction 
was  also  done  for  humidity  daily  averages.  The  analysis 
was  similar  to  that  in  figure  10,  except  that  the  Y-axis 
was  reversed. 


Sunrise  Sunset 
(Minimum  Tair)  (Maximum  Tair) 

(Maximum  RH)  (Minimum  RH) 


Figure  10.— A  diagram  analyzing  the  use  of  a  sine  function  to  ap- 
proximate average  air  temperature  and  relative  humidity  during 
daylight  hours,  given  only  maximum  and  minimum  values. 


RH     =  RH  -  0.212  (RH  -  RH  ). 

ave  *■  mm' 

The  regression  between  hourly  averaged  RH  and  the 
sine  function-weighted  daily  RH  was 

RHave  =  -17.0  +  1.32  RHsf   R2  =  0.84  n  =  120  days 

This  corrected  RH  was  also  an  improvement  over  a 
daily  RH  produced  merely  by  averaging  the  maximum 
and  minimum  RH. 

Again,  it  should  be  stated  that  these  data  reductions 
must  be  done  outside  of  the  main  DAYTRANS  program. 
Also,  this  sine  function  approximated  diurnal  humidity 
curve  would  not  be  appropriate  when  using  an  absolute 
measure  of  atmospheric  water  vapor  such  as  dewpoint. 

Soil  Temperature — Z(5) 

Because  the  diurnal  temperature  variation  at  20-cm 
depth  is  usually  less  than  1°  C,  any  midday  soil  tempera- 
ture^ measurement  is  an  adequate  daily  average  for  the 
purposes  of  this  model. 

Incoming  Shortwave  Radiation — Z(6) 

If  total  daily  radiation  is  collected,  radiation  divided 
by  day  length  in  minutes  will  provide  the  daily  average 
used  in  Z(6).  Day  length  for  a  flat  surface  is  computed  in 

Z(7). 

In  the  absence  of  any  radiation  data,  daily  potential 
radiation  can  be  determined  (Gamier  and  Ohmura  1968, 
Swift  1976).  Assuming  a  sine  wave  approximation  of  a 
cloudless  diurnal  radiation  trace,  peak  daily  radiation 
times  0.7  provides  a  fair  approximation  of  daily  average 
incoming  radiation  for  relatively  clear  areas. 

Day  Length— Z(7) 

XD     =  JD  -  79 

IF(XD  .  LT  .  0.0)  XD  =  286.0  +  JD 
DAY  =  3.125-(SIN(XD-0.01721))  +  12.0 
Z(7)     =  DAY  •  3600  ■  0.8 

This  sine  function,  driven  by  Julian  date  (JD),  com- 
putes day  length  in  hours  (DAY),  at  41°  N.  latitude,  for  a 
flat  surface,  on  any  day  of  the  year.  Because  this  day 
length  calculation  overestimates  the  duration  of  daylight 
with  sufficient  intensity  to  open  stomata  (0.1  ly  min-1), 
analysis  of  diurnal  radiation  traces  showed  that  a  0.8 
correction  more  closely  approximates  the  length  of  ef- 
fective daylight  for  transpiration.  The  sine  equation  can 
be  adjusted  for  any  latitude  by  the  3.125  coefficient  (3.5 
for  45°  N.)  and  is  typically  accurate  to  within  15  minutes 
per  day.  More  complex  equations  correcting  for  differ- 
ent slopes  and  aspects  are  also  available  (Gamier  and 
Ohmura  1968,  Swift  1976). 


AUXILIARY  CONSTANTS— B(i) 

Unit 

B(l)    =  maximum  soil  surface 

infiltration  rate  cm  day-1 

B(2)  =  root  zone  soil  capacity  cm3 
B(3)    =  tree  sapwood  storage  capacity  cm3 
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B(4) 

=  tree  or  stand  leaf  area 

cm2 

B(5) 

=  maximum  canopy  leaf 

conductance,  k, 

cm  s1 

B(6) 

=  canopy  and  litter  interception 

coefficient 

cm 

B(7) 

=  ground  surface  area 

cm2 

B(8) 

=  snowmelt  coefficient 

cm  °C-1 

B{9) 

=  midcrown  tree  height 

m 

B(10) 

=  spring  minimum  predawn  ^ 

MPa 

B(ll) 

=  stomatal  closure  threshold 

MPa 

INTERMEDIATE  VARIABLES— G(i) 

G(l)  =  effective  precipitation  (cm3  day-1) 

G(l)  -  AMAX1((Z(2)-(B(4)/B(7))  ■  B(6))  •  B(7),  0.0) 
where 
Z(2)  =  precipitation, 
B(4)  =  stand  leaf  area, 

B(6)  =  canopy  and  litter  interception  coefficient, 
B(7)  =  ground  surface  area. 

This  equation  reads  input  precipitation,  Z(2),  sub- 
tracts canopy  and  litter  interception,  B(6),  as  a  function 
of  leaf  area  index,  B(4)/B(7),  and  multiplies  by  ground 
surface  area,  B(7),  to  give  volume  of  water  input.  See 
G(50)  in  H20TRANS. 

G(2)  =  available  water  fraction,  root  zone  (cm3  cm-3) 
G(2)  =  X(2)/B(2) 

where 

X(2)  =  current  root  zone  water  storage, 
B(2)  =  root  zone  water  capacity. 

G(3)  =  available  water  fraction,  sapwood  (cm3  cm-3) 
G(3)  =  X(3)/B(3) 

where 

X(3)  =  current  sapwood  water  storage, 
B(3)  =  sapwood  storage  capacity. 

G(4)  =  snowmelt  (cm3  day-1) 

G(4)  =  AMAXl(B(8)-Z(3)-B(7),  0.0) 
IF  (X(l)  .LE.  0.0)  G(4)  =  0.0 

where 

B(8)  =  snowmelt  coefficient, 

Z(3)  =  air  temperature, 

B(7)  =  ground  surface  area, 

X(l)  =  current  snowpack  water  content. 

This  equation  melts  the  snowpack,  X(l)  at  a  rate  B(8), 
based  on  average  day  air  temperature,  Z(3),  multiplied 
by  ground  surface  area,  B(7),  and  gives  the  volume  of 
water  melted.  See  G(51)  in  H20TRANS. 

G(5)  =  leaf  area  index  (cm2  cm-2) 

G(5)  =  B(4)/B(7) 

where 
B(4)  =  stand  leaf  area, 
B(7)  =  ground  surface  area. 

G(9)  =  canopy  and  litter  evaporation  (cm3  day-1) 


G(9)  =  Z(2)-B(7)-G(l) 

where 

Z(2)  =  precipitation, 

B(7)  =  ground  surface  area, 

G(l)  =  effective  precipitation  into  the  soil. 

Evaporation  of  water  from  canopy  and  litter  intercep- 
tion is  calculated  by  G(9). 

G(10)  =  predawn  leaf  water  potential,  Bi^,  (-MPa) 

G(10)  =  AMAX1(B(10),  0.2/G(2)  +  0.01-B(9),  0.0) 
where 

B(10)  =  spring  minimum  predawn 

0.2     =  curve-fitting  coefficient, 

G(2)    =  available  water  fraction,  root  zone, 

0.01    =  hydrostatic  gradient  constant  (MPa  m_1), 

B(9)    =  midcrown  tree  height. 

Pre-dawn  leaf  water  potential  is  predicted  from  avail- 
able soil  moisture  G(2),  and  is  corrected  for  tree  height. 
See  G(10)  in  H20TRANS. 

G(12)  =  light-saturated  leaf  area  index  (m2  m-2) 
G(12)  =  10.0-(1.0-EXP(-4.6-Z(6))) 

where 

10.0,  1.0,  -4.6  =  curve-fitting  coefficients, 
Z(6)  =  incoming  shortwave  radiation. 

This  equation  computes  the  amount  of  leaf  area  index 
that  would  be  receiving  at  least  10%  of  full  sunlight  at 
the  measured  radiation  intensity,  Z(6),  and  assumed  ver- 
tical attenuation.  See  G(ll)  in  H20TRANS. 

G(13)  =  radiation-induced  k,  limiting  fraction  (dim.) 

G(13)  =  1.0-((G(5)-G(12))/G(5)) 
IF(G(12)  .GT.  G(5))  G(13)  =  1.0 

where 

G(5)    =  leaf  area  index, 

G(12)  =  light-saturated  leaf  area  index. 

Canopy  leaf  conductance  is  reduced  by  G(13),  the 
fraction  of  the  canopy  receiving  insufficient  radiation  to 
open  stomata.  See  G(12)  in  H20TRANS. 

G(14)  =  k,,  predawn  leaf  water  potential  correction 
(cm  s_1) 

G(14)  -  B(5)  -  (B(5)/B(11)-B(10))-(G(10)-B(10)) 
IF(G(10)  .GE.  B(ll))  G(14)  =  0.005 

where 

B(5)  =  maximum  canopy  k,, 

B(5)/B(ll)-B(10)  =  the  slope  of  the  k,  reduction  calcu- 
lated by  the  range  of  k,  (maximum 
(B(5))  to  0)  divided  by  the  corre- 
sponding range  of  ^  (spring  mini- 
mum \pv  B(10)  to  stomatal  closure 
point,  B(ll)j. 

G(10)  =  predawn  leaf  water  potential, 

B(10)  =  spring  minimum  \pv 

0.005  =  minimum  k,,  stomatal  closure. 

This  important  function  relates  reduction  in  canopy  k, 
to  soil  moisture  deficits  through  the  predawn  leaf  water 
potential  function,  G(10).  See  the  discussion  in 
H20TRANS  for  G(13). 
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G(15)  =  kp  humidity  correction  (cm  s_1) 


See  G(76)  in  H20TRANS. 


G(15)  =  G(14)-(G(14)-0.05-(Z(4)-106-4.0)) 


IF(G(15)  .LE.  0.0)  G(15)  =  0.005 


where 

G(14) 

=  kp  leaf  water  potential  effects, 

0.05 

=  slope  of  k,  reduction, 

Z(4) 

=  absolute  humidity  deficit, 

106 

=  unit  scalar, 

4.0 

=  humidity  control  threshold  factor, 

0.005 

=  stomatal  closure  constant. 

This  function  produces  the  humidity  effects  on  k, 
analagous  to  G(14)  in  H20TRANS. 

G(16)  -  kp  air  temperature  correction  (cm  s-1) 


G(16)  =  G(15)  +  G(15)- 0.003 -(Z(3)-10.0) 
IF(Z(8)  .LT.  0.0)  G(16)  =  AMAX1(G(15) 
+  0.02-Z(8),  0.005) 


where 

G(15) 

kp  humidity  effects, 

0.02,  0.003  = 

slope  of  k,  reduction  (cm  s_1  °C_1 

Z(3) 

air  temperature, 

10.0 

temperature  reduction  threshold, 

Z(8) 

night  minimum  air  temperature, 

0.005 

k  at  stomatal  closure. 

These  equations  increase  k,  slightly  with  air  tempera- 
tures above  10°  C,  reduce  k,  at  temperatures  between  0° 
and  10°,  and  markedly  reduce  k,  at  night  minimum  tem- 
peratures below  0.0°.  See  special  function  Si, 
H20TRANS. 

G(17)  =  kp  final  canopy  average  (cm  s^1) 
G(17)  =  G(16)-G(13) 

where 

G(16)  =  kp  plus  humidity,  air  temperature  and  leaf 

water  potential  controls, 
G(13)  =  kp  radiation  reduction  fraction, 
See  G(17)  in  H20TRANS. 

G(18)  =  transpiration  flux  density  (g  cm-2  s_1) 
G(18)  =  Z(4)-G(17) 

where 

Z(4)    =  absolute  humidity  deficit, 
G(17)  =  kp  final  canopy  average, 
See  G(17)  in  H20TRANS. 

G(19)  =  daily  stand  transpiration  (g  day1) 
G(19)  =  G(18)-B(4)-Z(7) 

where 

G(18)  =  transpiration  flux  density, 

B(4)    =  stand  leaf  area, 

Z(7)    =  day  length, 

Final  stand  transpiration  G(19)  =  F(3,4). 

G(20)  =  transpiration/potential  transpiration  ratio  (dim.) 
G(20  =  G(17)/B(5) 

where 

G(17)  =  kp  final  canopy  average, 
B(5)    =  maximum  kr 


G(21)  =  areal  transpiration  (mm  day1) 
G(21)  =  G(19)/B(7)-10.0 

where 

G(19)  =  daily  stand  transpiration, 
B(7)    =  ground  surface  area. 

G(50)  =  precipitation  and  snowmelt  water  input  (cm3 
day1) 

G(50)  =  G(4) 
IF  (X(l)-G(4)  .LT.  0.0)  G(50)  =  X(l) 

where 

G(4)    =  snowmelt, 

X(l)    =  snowpack  water  content. 

Water  input  into  X(2)  comes  from  precipitation  and 
snowmelt.  Snowmelt,  G(4),  cannot  exceed  the  current 
snowpack  water  content,  X(l).  G(50)  equals  F(l,2). 

G(51)  =  surface  runoff  (cm3  day-1) 
G(51)  =  0.0 

IF  ((G(l)  +  G(4))/B(7)  .GT.  B(l))  G(51)  =  (G(50)-B(l))-B(7) 


where 

G(l) 

=  precipitation, 

G(4) 

=  snowmelt, 

B(7) 

=  ground  surface  area, 

B(l) 

=  maximum  soil  surface  infiltration  rate, 

G(50) 

=  precipitation  and  snowmelt  water  input 

Water  input  in  excess  of  the  maximum  soil  surface  in- 
filtration rate  is  routed  to  surface  runoff.  G(51)  equals 
F(2,5). 

G(52)  =  subsurface  outflow  (cm3  day"1) 
G(52)  =  0.0 

IF  (X(2)  +  G(50)  .GT.  B(2))  G(52)  =  X(2)  +  G(50)-G(51)-B(2) 
where 

X(2)    =  current  root  zone  soil  water  content 

G(50)  =  precipitation  and  snowmelt  water  input, 

G(51)  =  surface  runoff, 

B(2)    =  root  zone  soil  water  capacity. 

If  water  input,  G(50),  minus  surface  runoff  exceeds 
the  capacity  of  X(2),  the  excess  water  goes  to  subsurface 
outflow.  G(52)  equals  F(2,6). 

G(53)  =  root  water  uptake  (cm3  day-1) 

G(53)  =  B(3)-X(3)  +  G(19) 

where 

B(3)    =  tree  sapwood  storage  capacity, 
X(3)    =  current  sapwood  water  storage, 
G(19)  =  daily  transpiration. 

This  equation  computes  water  transfer  from  X(2)  to 
X(3)  based  solely  on  deficits  in  X(3)  and  transpiration  de- 
mand. In  a  previous  model  version,  I  had  a  soil-water  up- 
take resistance  also  controlling  water  movement  at  this 
point,  but  it  seems  to  be  unnecessary. 

The  feedback  resistance  produced  by  the  effects  of 
soil  water  deficit  on  predawn  leaf  water  potential, 
G(10),  leaf  conductance,  G(13),  and  transpiration,  G(19), 
controls  the  overwithdrawal  of  water  from  X(2)  except 
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after  extreme,  extended  water  stress.  At  that  point  the 
model  should  predict  tree  mortality! 

Two  additional  capabilities  were  recently  added  to 
DAYTRANS,  although  they  were  not  included  in  the 
results  section. 

1.  Penman-Monteith  calculation  of  transpiration 

FUNCTION  PENMON(TAIR,RAD,VPD,XLAI,COND) 
GAMMA  =  0.646  +  0.0006*TAIR 
PLAI  =  XLAI/2. 
Tl  =  TAIR  +  0.5 
T2  =  TAIR-0.5 

SVP1  -  6.1078*EXP((17.269*Tl)/(237.0  +  Tl)) 
SVP2  =  6.1078*EXP((17.269*T2)/(237.0  +  T2)) 
SLOPE  =  SVP1-SVP2 
XNETR  =  RAD*0.8*697.3 
CP=1.0lE  +  3 
PA  =  1.292-0.00428*TAIR 
RAZO- 
RS =  100./COND 

XLAT  =  (2.501  -  0.0024*TAIR)  *  1.0E  +  6 
XTRANS  =  ((SLOPE*XNETR)/PLAI 
C+(CP*PA)*(VPD/RA))/ 
C  (SLOPE  +  GAMMA*(1.0  +  RS/RA)) 
PENMON  =  XTRANS/(XLAT  *  10.) 
RETURN 
END 

This  is  an  adaptation  of  the  widely  used  Penman- 
Monteith  equation  (Jarvis  et  al.  1976): 


(AR17/PLAI)  +  Cp  P(VPD)/ra 
XE  =  — 


A  +  7  (1  +  r/r) 


where 
X 

E 

PENMON 
A 

Rrj 

PLAI 

VPD 
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XLAT  =  latent  heat  of  vaporization  of 
water  (J  kg-1), 

XTRANS  =  evaporation  (Wm-2), 
evaporation  (g  cm-2  s_1), 
SLOPE  =  slope  of  the  saturation  vapor 
pressure  curve  at  Tair,  (mbar  °C_1), 
XNETR  =  net  radiation  (W  m"2), 
projected  leaf  area  index, 
CP  =  specific  heat  of  air  (J  kg ~ 1  °C_1), 
PA  =  density  of  air  (Kg  m~3), 
vapor  pressure  deficit  (mb), 
RA  =  aerodynamic  resistance  (s  m"1), 
GAMMA   =    psychrometric  constant 
(mb  °C-1),  calculated  by  7  =  CpP/X, 
rs  =  RS  =  stomatal  resistance,  calculated  as 

1/k,  G(17),  (s  m-1)). 

Under  conditions  of  high  net  radiation  but  low  air 
temperature  common  in  spring,  radiant  energy  becomes 
a  significant  driving  variable  for  evaporation.  This 
equation  assumes  a  constant  ra  of  5.0  s  m_1  with  no 
windspeed  input.  It  also  divides  the  net  radiation  term 
by  projected  leaf  area  index  to  more  realistically  ex- 
press the  radiant  energy  loading  on  different  canopy 
layers.  This  calculation  of  stand  transpiration  can  be 
substituted  for  G(17)  in  H20TRANS  or  G(18) 
DAYTRANS. 


2.  Net  photosynthesis  subroutine 
SUBROUTINE  PHOTO 
COMMON  K,X(2,20),F(20,20),G(100), 
B(100),Z(20),Y(160) 

C 

C  Z(10)  =  CANOPY  AVERAGE  RADIATION 

C 

Z(10)  =  (Z(6)  +  Z(6)*EXP(-0.7*G(5)/2.))/2. 

C 

C  Z(9)  =  AVERAGE  NIGHT  TEMPERATURE 
C 

Z(9)  =  (Z(3)  +  Z(8))/2. 

C 

G(64)  =  ((Z(10)-0.0143)/(Z(10)  +  0.322))* 
C  (0.0182  +  0.0105*Z(3)-0.000194*(Z(3)* *2)) 

IF(G(64).LT.0.)G(64)  =  0. 

G(65)  =  (0.0006*(G(17)/1.6)*G(64))/ 
C  ((G(17)/l.6)  +  G(64)) 

G(62)  =  0.001  *(24.-Z(7)/3600.)*EXP(0.2*Z(9)) 

G(66)  =  G(65)*Z(7) 

G(67)  =  G(66)-G(62) 

RETURN 

END 

This  subroutine  was  derived  from  the  photosynthesis 
model  of  Lohammar  et  al.  (1980). 

G(65)  =  PSN  =  gross  photosynthesis  (mg  cm-2  s_1) 

AC02(k,-gJ 


G(65)  =  PSN  = 


where 

AC02  =  COz  gradient  from  the  atmosphere  to  the 
carboxylation  site  in  the  leaf  (=0),  taken  as 
0.0006  kg  m-3 

k,  =  G(17)  =  leaf  conductance  to  water  vapor 
(cm  s_1), 

1.6      =  ratio  of  the  diffusion  coefficients  of  water 

vapor  and  CO,  in  the  air, 
gm       =  G(64)  =   mesophyll  C02  conductance 

(cm  s- '). 

G(64)  is  calculated  by  Lohammar  et  al.  (1980)  as: 
G(64)  =  gm  =-I-4L-cxf(T) 


where 
I 


am 


I  +  h 

Z(10)    =    incoming    shortwave  radiation 

(W  m~2)  with  a  conversion  factor  of  1  W  m~2 

=  1.43  x  10  3  cal  cm-2mh-r\ 

light  compensation  point,  where  net  PSN  = 

0,  given  as  10  W  m"2  =  0.0143  cal  cm  2 

min~\ 

the  irradiance  at  which  g  is  half  of  its  max- 
imum  value,  given  as  225  W  m  2  =  -0.322 
cal  cm-2  min- \ 

the  maximum  gm  at  any  air  temperature  T, 
generated  by  the  polynomial  equation: 


in 


a  =  0.0182  +  0.0105  (Z(3))  +  0.000194  (Z(3)2) 

This  curve  was  fitted  to  data  from  Nielson  and  Jarvis 
(1976),  R2  =  0.67,  and  checked  against  a  variety  of  other 
sources  for  conifers.  It  calculates  a  maximum  gm  =  0.16 
cm  s"1  at  26.3°  C. 


17 


G(66)  =  Daily  net  C02  uptake  (mg  cm"2) 
G(66)  =  G(65)  •  Z(7) 

where 

G(65)  =  gross  photosynthesis  (mg  cm~2  s1), 
Z(7)    =  daylength  (s). 

G(62)  =  Night  foliar  respiration  loss  (mg  cmr2 

G(62)  =  0.001  •  (24  -  (Z(7)/3600))  ■  exp  (0.2  •  Z(9)) 
where 

24  -  (Z(7)/3600)  =  night  length  (h), 
Z(9]  =  average    night    air  temperature 

(°C). 

This  equation  was  used  by  Emmingham  and  Waring 
(1977). 

Finally,  24-hour  net  photosynthesis  is  approximated 
by 

G(67)  =  24-hour  net  C02  uptake  (mg  cm  2) 
G(67)  =  G(66)  -  G(62) 

where 

G(66)  =  daily  net  C02  uptake, 
G(62)  =  night  foliar  respiration  loss. 


RESULTS 

This  preliminary  validation  had  two  objectives.  First 
was  to  test  if  H20TRANS  and  DAYTRANS  were  capable 
of  simulating  the  diurnal  and  seasonal  patterns  of 
transpiration  and  water  stress  development  in  lodgepole 
pine.  Because  transpiration  of  a  stand  cannot  be  meas- 
ured directly,  validation  focused  on  diurnal  and 
seasonal  patterns  of  measurable  variables  ^  and  k,  in 
individual  trees  and  soil  moisture  depletion.  The  second 
objective  was  to  determine  if  the  much  smaller,  simpli- 
fied DAYTRANS  model  could  predict  seasonal  patterns 
with  minimal  resolution  loss  compared  to  H20TRANS. 
The  H20TRANS  model  was  designed  as  a  research  tool 
and  requires  more  detailed  climatic  data  and  special- 
ized parameter  estimates.  DAYTRANS  was  an  attempt 
to  scale  the  logic  of  the  H20TRANS  model  down  to  a 
level  where  input  data  requirements  would  be  feasible 
for  operational  use. 

FIELD  VALIDATION  SITE 

During  the  fall  of  1977,  a  study  site  was  established  at 
the  USDA  Forest  Service  Fraser  Experimental  Forest  in 
the  central  Colorado  Rocky  Mountains.  The  site  studied 
supported  an  uneven-aged  Pinus  contorta  stand,  with  oc- 
casional Populus  tremuloides,  ranging  in  age  from  10 
years  to  the  12-m  tall  canopy  dominants  at  60  years.  The 
elevation  was  2,700  m,  and  topography  level  on  a  glacial 
outwash.  Stand  density  was  2,740  trees  per  hectare 
with  a  basal  area  of  30.7  m2  ha-1  and  site  index  of 
22.6  m  at  100  years. 

The  climate  of  this  area  is  cool  and  dry  for  coniferous 
forests.  Temperature  extremes  of  from  -40°  C  to  32°  C 
have  been  recorded.  Frost  is  possible  on  any  night  dur- 


ing the  growing  season.  Minimum  temperatures  on  site 
during  the  summer  of  1978  were  always  below  4°  C,  and 
maximum  daily  air  temperatures  averaged  from  20°  C  to 
23°  C,  with  a  high  of  26°  C.  Midday  relative  humidity 
typically  ranged  from  15%  to  25%,  except  during  thun- 
derstorms. Precipitation  averages  58.4  cm  per  year, 
with  nearly  two-thirds  falling  as  snow  between  October 
and  May  (Alexander  and  Watkins  1977). 

Beginning  on  Julian  date  130  (May  10),  1978,  con- 
tinuous meteorological  data  was  recorded  until  JD  250 
(September  10),  1978.  Air  temperature  measured  by  a 
thermistor  was  continuously  recorded  by  strip  chart. 
Dewpoint  temperature  was  measured  with  a  heated  lith- 
ium chloride  dew  point  sensor,  also  recorded  continu- 
ously on  the  strip  chart.  These  instruments  were 
mounted  in  a  standard,  vented  weather  station  box 
placed  1.5  m  above  the  ground.  Onsite  soil  temperature 
at  15-cm  depth  was  taken  with  a  dial  temperature  probe 
once  a  day.  Also,  a  sling  psychrometer  was  used  every 
two  hours  during  daylight  hours  to  take  wet-  and  dry- 
bulb  temperatures  for  backup  measurement  of  air  tem- 
perature and  humidity.  Precipitation  and  a  continuous 
trace  of  incoming  shortwave  radiation  were  being 
recorded  5  km  away  at  the  Fraser  Experimental  Forest 
headquarters.  At  the  end  of  the  summer,  these  data 
were  compiled  into  the  format  required  for  input  into 
H20TRANS  and  DAYTRANS  as  Z,  or  driving,  variables. 

PARAMETER  ESTIMATES 

To  run  the  models,  a  number  of  site  and  stand 
parameters  were  also  needed.  Because  a  complete  bio- 
mass  survey  of  the  stand  could  not  be  done  and  valida- 
tion was  focusing  on  individual  tree  responses,  a  single 
tree  representative  of  the  stand  was  defined  for  model- 
ing. The  characteristics  of  this  tree  follow: 

Ground  surface  area  =  3.5  m2 
Diameter  (b.h.)  =  8.2  cm 
Height  =  7.0  m 
Needle  weight  =  2.87  kg 
Needle  area  =  2.1  x  105  cm2 
LAI  =  6.0 

Sapwood  basal  area  =  73.9  cm2 
Sapwood  volume  =  21,000  cm3. 
Trees  of  this  general  size  had  been  used  to  gather  data 
for  some  of  the  physiological  responses  in  the  model. 
Consequently,  these  parameter  estimates  can  be 
verified  from  measured  biomass  data  (Running  1980c). 
Also,  the  "observed"  field  data  used  in  this  preliminary 
validation  were  taken  on  trees  of  this  size. 

A  list  of  initial  conditions  for  the  X(i),  state  variables, 
at  the  beginning  of  the  H20TRANS  and  DAYTRANS 
model  runs  on  JD  130,  1978,  and  the  B(i),  constant 
parameters,  is  given  in  table  1. 

VALIDATION  RUNS 

Simulation  runs  were  made  for  JD  130  to  245,  1978,  by 
H20TRANS  and  DAYTRANS.  The  overall  water  budgets 
predicted  by  the  two  models  were: 
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H20TRANS  DAYTRANS 


Soil  depletion 


Transpiration 


Transpiration 
Evaporation 
Subsurface  outflow 
Total 


36% 
15% 
49% 
42.9  cm 


38% 
17% 
45% 


This  was  taken  from  a  total  water  supply  during  the 
period  consisting  of  snowpack,  soil  water,  and  precipita- 
tion. Figure  11  plots  the  seasonal  progression  of 
transpiration  and  soil  moisture  depletion  simulated  by 
the  models,  and  the  actual  soil  moisture  depletion  meas- 
ured on  site  by  technicians  using  a  neutron  probe.  The 
neutron  probe  data,  taken  at  monthly  intervals,  from  0 
to  1.5  m  depth,  at  15-cm  increments,  from  four  access 
tubes,  is  shown  in  Running  (1980a).  The  initiation  of  soil 
moisture  depletion  and  final  magnitude  of  drawdown 
were  predicted  fairly  accurately  by  both  models.  The 
models  appear  to  overestimate  water  depletion  rates 
somewhat  during  midsummer;  however,  the  slope  of 
overall  seasonal  depletion  was  fairly  closely 
represented.  DAYTRANS  predicted  transpiration  rates 
roughly  5%  higher  than  H20TRANS  during  midsummer 
beginning  around  JD  170.  This  discrepancy  is  attributed 
to  the  difficulties  of  averaging  evaporative  demand  ac- 
curately for  the  day  and  scaling  the  humidity  reduction 
of  k,  correctly  in  DAYTRANS  to  averaged  values.  Final 
seasonal  transpiration  estimates  agree  to  within  4%, 
which  is  remarkable  considering  the  much  reduced 
resolution  of  DAYTRANS. 


■DAYTRANS 
-H2OTRANS 
-  Observed 


130 


Julian  date 


Figure  11.— Simulation  runs  comparing  H20TRANS  and  DAYTRANS 
results  to  observed  data  on  lodgepole  pine  at  the  Fraser  Ex- 
perimental Forest  study  site  during  the  summer  of  1978.  Soil 
water  depletion  was  measured  by  a  neutron  probe.  There  is  no  ob- 
served data  for  transpiration;  therefore,  only  the  comparison  of 
model  predictions  is  shown. 

Predawn  leaf  water  potential,  B^,,  is  considered  by 
many  tree  physiologists  to  be  the  best  and  most  easily 
measured  way  of  assessing  seasonal  water  stress 
development  in  a  tree.  In  figure  12,  the  H20TRANS  and 
DAYTRANS  predictions  of  seasonal  stress  development 
BiA,  are  shown  compared  to  a  compilation  of  observed 
B^,  measurements  taken  at  different  lodgepole  pine 
trees  at  the  Fraser  site  throughout  the  summer  of  1978. 
Again,  both  models  performed  well  relative  to  each 
other  and  the  measured  field  data. 


Table  1.— Model  conditions  for  the  preliminary  validation  exercise 


H20TRANS 


DAYTRANS 


X(i):  Initial  conditions  (cm3) 


X(1) 

=  7.2  X 

105 

X(2) 

=  0.0 

X(3) 

=  1.0  x 

104 

X(4) 

=  3.0  x 

105 

X(5) 

=  3.0  x 

105 

X(6) 

=  0.0 

X(7) 

=  1.6  x 

103 

X(8) 

=  0.9  x 

102 

X(9) 

=  0.0 

X(10) 

=  0.0 

B(i)  Constant  parameters 

B(1) 

=  0.015 

cmh"' 

B(2) 

=  10.0 

cm  h~ 

B(3) 

=  3.5  x 

104 

cm3 

B(4) 

=  3.0  x 

105 

cm3 

B(5) 

=  3.0  x 

105 

cm3 

B(6) 

=  0.010 

cm  LAI 

B(7) 

=  1.6  x 

104 

cm3 

B(8) 

=  9.0  x 

102 

cm3 

B(9) 

=  2.1  x 

105 

cm2 

B(10) 

=  3.5  x 

104 

cm2 

B(11) 

=  7.0 

m 

B(12) 

=  0.15 

cm  s- 1 

B(13) 

=  0.4 

MPa 

B(14) 

=  1.5 

MPa 

B(15) 

=  74.0 

cm2 

B(16) 

=  2.0 

h"1 

B(17) 

=  20.0 

h-1 

X(i):  Initial  conditions  (cm3) 

X(1)  =  7.2  x  105 

X(2)  =  6.0  x  105 

X(3)  =  1.33  x  104 

X(4)  =  0.0 

X(5)  =  0.0 

X(6)  =  0.0 

X(7)  =  0.0 


B(i):  Constant  parameters 


B(D 

5.0  x 

101 

cm  day" 

B(2) 

6.0  x 

105 

cm3 

B(3) 

1.33  x  104 

cm3 

B(4) 

2.1  x 

105 

cm2 

B(5) 

0.15 

cm  s_1 

B(6) 

0.30 

cm  LAI" 

B(7) 

3.5  x 

104 

cm2 

B(8) 

0.15 

cm  °C~ 

B(9) 

7.0 

m 

B(10) 

0.4 

MPa 

B(11) 

1.5 

MPa 
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Leaf  conductance  is  the  most  easily  measured  physio- 
logical variable  that  is  closely  related  to  transpiration. 
To  evaluate  the  model's  predictions  of  k,,  the  seasonal 
progression  of  maximum  daily  k,  was  plotted  in  figure  13 
and  compared  to  observed  k,  data  compiled  from 
measurements  on  different  lodgepole  pine  taken 
throughout  the  summer  of  1978.  The  beginning  of  stress 
development  around  JD  185  is  evident  in  figure  12,  with 
corresponding  reducing  of  k,  in  figure  13.  Note  that  in 
both  figures  12  and  13  the  models  are  relatively  success- 
ful in  predicting  (1)  the  time  of  initial  stress  develop- 
ment, and  (2)  the  overall  rate  of  increase  in  stress  (fig. 
12)  and  reduction  in  k,  (fig.  13). 

To  further  examine  the  ability  of  H20TRANS  to  simu- 
late diurnal  dynamics  of  \px  and  k,,  19  days  of  ^  and  k, 
data  from  12  different  trees  at  the  Fraser  site  were  com- 
piled. Each  data  point,  representing  a  specific  date  and 
time  of  day,  was  plotted  against  the  corresponding 
H20TRANS  simulated  value  for  that  date  and  time.  The 
results  of  this  analysis,  shown  in  figure  14  for  i/',  and 
figure  15  for  k,,  cover  an  extensive  array  of  conditions  of 
temperature,  humidity,  radiation,  and  soil  moisture  ex- 
pected in  measurements  taken  from  sunrise  to  sunset 
throughout  the  summer.  Both  the  slope  and  Y-intercept 
of  the  linear  regression  equation  in  figure  14  correspond 
well  with  the  optimal  1:1  where  every  measured  value  is 
predicted  exactly.  It  is  evident  in  analyzing  the  scatter 
of  points  that  the  model  tends  to  underestimate  ^  in  the 
1.0-1.3  MPa  range  and  overestimate  ^  in  the  1.3-1.5 
MPa  range. 

The  diurnal  prediction  of  k,  was  less  successful,  possi- 
bly because  of  the  large  number  of  factors  controlling  k. 
The  Y-intercept  of  the  regression  line  shows  that  at  low  k, 
values  (below  0.06  cm  s1)  H20TRANS  often  over- 
estimates k,.  At  k,  values  above  0.08  cm  s  1  the  model 
underestimates  k,.  This  illustrates  the  concept  that  when 
average  "conservative"  response  equations  are  used  in  a 
model,  the  deviations  from  the  average  will  be  most 
noticeable  at  the  extremes  of  the  response  curves.  For- 
tunately, these  high  resolution  errors  compensated  each 
other,  allowing  reasonable  estimates  of  the  longer  term 
seasonal  trends  of  tree  water  use  and  stress  development. 
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Figure  12.— A  comparison  of  H20TRANS  and  DAYTRANS  predic- 
tions of  increase  in  tree  water  stress  (predawn  leaf  water  poten- 
tial, Bi/,)  with  measured  values  on  lodgepole  pine  at  the  Fraser 
Experimental  Forest  study  site  during  the  summer  of  1978. 
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Figure  13.— Predicted  and  observed  levels  of  maximum  daily  leaf 
conductance,  k,,  on  lodgepole  pine  at  the  Fraser  site  throughout 
the  summer  of  1978.  Observed  points  represent  an  average  of 
9-11  diffusion  porometer  measurements  taken  throughout  the 
tree  crown. 

SENSITIVITY  ANALYSIS 

To  assess  the  response  of  the  overall  model  to  various 
"perturbations,"  a  sensitivity  analysis  was  run  on 
DAYTRANS.  A  similar  test  was  not  run  on  H20TRANS 
because  the  models  are  so  similar  and  the  cost  and  com- 
puter time  of  running  the  larger  model  was  prohibitive. 
To  execute  this  analysis,  six  independent  variables 
were  individually  changed;  first  an  increase  of  50%, 
then  a  decrease  of  50%,  for  a  total  of  12  model  runs.  The 
variables  changed  were  initial  snowpack  water,  X(l), 
rooting  zone  water  capacity,  B(2),  leaf  area,  B(4),  max- 
imum k,,  B(5),  interception  coefficient,  B(6),  and  snow- 
melt  coefficient,  B(8).  Three  dependent  variables  were 
analyzed  to  determine  changes  brought  about  by  these 
+  /-50%  perturbations.  The  dependent  variables  were 
total  transpiration,  X(4),  total  subsurface  outflow,  X(6), 
and  predawn  leaf  water  potential,  G(10).  The  results  of 
this  sensitivity  test  are  given  in  table  2. 

Changes  in  the  snowpack  primarily  caused  changes  in 
subsurface  outflow.  Snowmelt  goes  through  the  system 
relatively  quickly  in  the  spring;  therefore,  once  the 
rooting  zone  soil  is  saturated,  the  water  is  gone  before 
trees  can  use  much.  The  rooting  zone  is  the  fundamental 
source  of  water  for  the  tree.  Consequently,  changes  in 
root  zone  water  capacity  had  large  effects  on  transpira- 
tion and  tree  water  stress  development.  This  is  probably 
the  most  difficult  parameter  in  the  model  to  accurately 
define.  Increasing  leaf  area  made  little  difference  on 
anything  except  causing  higher  water  stress  develop- 
ment, because  the  soil  water  was  depleted  sooner. 
Decreased  leaf  area  markedly  decreased  transpiration 
and  water  stress  and  increased  subsurface  outflow 
somewhat. 

Increased  maximum  k,  made  surprisingly  little  differ- 
ence to  transpiration;  water  was  used  more  quickly,  but 
the  overall  available  amount  changed  very  little.  This 
analysis  does  not  reflect  results  that  would  be  obtained 
from  different  sites.  On  a  less  water  stressed  site  that 
retained  substantial  available  soil  water  at  the  end  of 
the  summer,  increased  maximum  k,  or  increased  leaf 
area  would  cause  much  higher  transpiration. 
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Table  2.— A  sensitivity  analysis  of  DAYTRANS 


Dependent  variables 


Independent  variables 


X(4) 
Transpiration 


X(6)  G(10) 
Subsurface  Predawn 
outflow        leaf  y 


X(1):  Snowpack  water 
+  50% 
-  50% 

B(2):  Root  zone  water 
+  50% 
-50% 


B(4):  Leaf  area 


+  50% 
-  50% 


B(5):  Maximum  k. 

+  50% 
-50% 


B(6):  Interception  coeff. 

+  50% 
-  50% 

B(8):  Snowmelt  coeff. 

+  50% 
-50% 


+  2% 
-4% 


+  22% 
-  45% 


+  6% 
-36% 


+ 10% 
-43% 


-2% 
+  7% 


+  50% 
-  49% 


2% 
•7% 


-6% 
+  3% 


0 

+  5% 


+  2% 
-7% 


-3% 
+  4% 


-  48% 
+  40% 


+  25% 
-63% 


+  21  % 
-  66% 


+  2% 
- 10% 


The  results  from  the  change  in  interception  coeffi- 
cients are  rather  artificial  because  virtually  no  rain  fell 
during  the  summer  of  1978.  During  preliminary  model 
runs,  it  was  found  that  changes  in  interception  can 
radically  change  transpiration  and  water  stress  predic- 
tions. A  very  low  interception  rate  can  cause  a  small 
shower  to  recharge  the  soil,  eliminating  all  water  stress, 
while,  in  fact,  it  requires  a  rainfall  event  on  the  order  of 
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Figure  14.— A  summary  of  all  leaf  water  potentials,  y,,  measured 
with  a  pressure  chamber  on  lodgepole  pine  at  the  Fraser  site  dur- 
ing 1978  compared  to  predicted  ^  from  the  H20TRANS  simula- 
tion run.  The  points  represent  data  taken  at  all  times  of  the  day 
from  sunrise  to  sunset.  The  dashed  line  shows  the  1:1  ratio. 
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Figure  15.— A  summary  of  all  measured  crown  leaf  conductance 
averages  from  lodgepole  pine  on  the  Fraser  site,  compared  to 
predicted  k,  from  the  1978  H20TRANS  simulation  run.  Each  meas- 
ured point  is  the  average  of  9-11  diffusion  porometer  readings 
taken  throughout  the  tree  crown  on  foliage  ranging  from  0  to  3 
years  old.  Sampling  was  done  periodically  from  sunrise  to  sunset. 

1.5  cm  in  24  hours  to  cause  substantial  soil  recharge  in 
the  Fraser  stand. 

The  sensitivity  of  the  snowmelt  coefficient  also  ap- 
pears to  be  low,  but  this  could  be  misleading.  A  site  with 
a  larger  snowpack  would  show  higher  sensitivity  to  the 
influence  of  the  melt  rate  in  partitioning  snowmelt  to 
transpiration  versus  subsurface  outflow. 


CONCLUSIONS 

The  preliminary  validation  of  H20TRANS  and 
DAYTRANS  demonstrates  that  the  models  are  both 
capable  of  estimating  reasonable  water  balances  for 
lodgepole  pine  (fig.  11).  The  models  also  appear  capable 
of  estimating  the  physiology  of  seasonal  water  stress 
development  in  lodgepole  pine  (figs.  12  and  13).  In  par- 
ticular, H20TRANS  is  able  to  track  measured  \pi  and  k, 
diurnally  over  a  wide  variety  of  conditions  (figs.  14  and 
15).  Much  of  the  models'  emphasis  is  on  the  control  of  k, 
by  temperature,  radiation,  humidity  and  soil  water 
supply. 

The  use  of  morning  maximum  leaf  conductance  as  a 
starting  point  for  diurnal  k,  calculations  seems  valuable. 
The  linkage  of  morning  maximum  k,  to  soil  water  status 
through  predawn  leaf  water  potential  gave  fairly  ac- 
curate results. 

The  calculation  of  i/',  in  H20TRANS  is  unique  and  its 
success  (fig.  15)  can  be  attributed  to  two  things.  First, 
the  flow  resistance  used  to  calculate  i/-,  has  both  a  soil 
temperature  and  soil  moisture  component.  Second,  the 
flow  that  generates  ^\  is  the  root  water  uptake,  not  the 
transpiration  rate  (which  other  models  have  used).  This 
model  configuration  allows  for  capacitance  in  the  tree 
hydraulic  system  and  yields  the  i/-,  diurnal  hysteresis 
that  has  been  observed  in  the  field  by  many  researchers. 
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The  generality  of  this  i/-,  and  k,  paradigm  remains  uncer- 
tain until  it  can  be  tested  on  other  species  under  dif- 
ferent conditions. 

A  primary  conclusion  from  this  work  has  been  that 
prediction  of  tree  water  stress  development  requires 
more  than  just  accurate  physiological  response  func- 
tions. Of  equal  importance  is  sound  representation  of 
hydrologic  inputs  and  storages.  The  system  inputs  and 
storages  control  the  timing  of  when  water  stress  begins 
in  the  trees.  The  physiological  processes  primarily  con- 
trol the  subsequent  rate  of  tree  water  stress  develop- 
ment. Interception  and  surface  evaporation  losses  must 
be  accurate  to  determine  effective  precipitation  input. 
Potential  errors  from  these  sources  were  minimized  in 
this  validation  because  the  Fraser  site  had  only  7.4  cm 
of  precipitation  during  the  entire  summer  of  1978.  On 
wetter  sites,  these  processes  could  profoundly  influence 
the  success  of  a  simulation.  Of  equal  importance  in 
model  performance  is  defining  the  volume  of  the  soil 
rooting  zone,  the  basic  water  storage  reservoir  for  the 
model.  This  requires  knowledge  of  tree  rooting  depth 
and  density,  averaged  for  the  trees  to  be  modeled. 

One  avenue  of  future  research  could  use  the  models  to 
estimate  some  of  these  more  difficult  parameters.  It  may 
be  easier  to  measure  transpiration  throughout  a  season, 
and  with  the  model  calculate  back  to  the  size  of  rooting 
zone  that  would  have  been  required  to  sustain  the 
observed  activity.  Similarly,  interception  and  evapora- 
tion losses  could  be  analyzed  indirectly  by  measuring 
the  duration  and  intensity  of  precipitation  that  is 
necessary  before  tree  water  stress  recovery  can  be 
detected. 

Another  critical  component  of  water  stress  modeling 
is  dependable  measures  of  biomass,  particularly  leaf 
area,  and  to  a  lesser  extent,  sapwood  storage  volume. 
Meteorological  driving  variables  are  critical  to  model 
performance.  It  became  clear  making  simultaneous  runs 
of  H20TRANSand  DAYTRANS  that  many  problems  can 
develop  when  averaging  meteorological  data  beyond  ac- 
ceptable limits.  Whereas  a  single  daily  average  soil 
temperature  is  fairly  reasonable,  a  daily  average  incom- 
ing radiation  value  is  less  so,  unless  generated  by  an  in- 
tegrating data-logger.  Even  then  the  opportunity  to 
model  diurnal  radiation  thresholds  and  day  length  is 
lost. 

With  any  modeling  effort,  decisions  of  temporal  and 
spatial  resolution  must  be  balanced  against  the  level  of 
accuracy  required.  If  high  resolution  predictions  of  diur- 
nal t/^and  k,  are  required,  H20TRANS  should  be  run  on 
no  larger  site  than  a  few  hectares  of  even-aged  trees,  of 
a  single  species,  on  homogeneous  topography.  For  more 
general  seasonal  estimates  of  transpiration  and  overall 
stress  development,  DAYTRANS  could  provide  approx- 
imate answers  over  many  hundreds  of  hectares  oc- 
cupied by  a  number  of  tree  species. 

Given  the  above  considerations,  these  models  have  a 
variety  of  potential  applications.  They  can  provide  a 
physiologically  sensitive  estimate  of  evapotranspiration 
(ET)  for  hydrologic  and  meteorologic  studies.  In  studies 
researching  processes  sensitive  to  tree  water  stress, 
these  models  may  provide  a  prediction  of  water  stress. 


For  example,  throughout  the  western  United  States 
forest  stands  periodically  endure  water  stress  sufficient 
to  severely  impede  photosynthesis.  Good  correlations 
between  model-predicted  site  water  stress  and  site  pro- 
ductivity (Grier  and  Running  1977,  Running  1981)  are 
likely.  There  may  be  a  significant  correlation  between 
certain  insect  and  disease  epidemics  and  stress  develop- 
ment in  the  host  trees  (Mattson  and  Addy  1975,  Waring 
and  Pitman  1980).  A  major  component  of  fire  danger 
rating  systems  is  moisture  content  of  the  living  fuels. 
Prediction  of  plant  moisture  content  must  entail  a  plant 
water  stress  model.  Quantitative  estimates  of  forest 
nutrient  cycles  require  knowledge  of  soil  and  tree  water 
movement. 

These  models  should  be  tested  on  forest  stands  of 
varying  species,  age  and  structure  under  different 
physical  and  climatic  conditions  before  any  general  use 
is  attempted.  A  more  extensive  validation  of  H20TRANS 
is  in  progress  on  eight  different  stands  in  southern 
Wyoming  (Knight  et  al.  1984). 
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APPENDIX  1.— Variable  Notations. 


B^,  predawn  leaf  water  potential,  MPa,  treated  as  positive  values  in 

these  models 
^,  leaf  water  potential,  MPa 

k,  leaf  conductance  to  water  vapor  diffusion,  cm  s  ~ 1 

LA  tree  or  stand  leaf  area,  cm2 

LAI  leaf  area  index,  the  ratio  of  canopy  leaf  area  to  ground  area 

ABSHD         absolute  humidity  deficit,  g  cm-3 

T  transpiration  rate,  cm3  h" 1 

TFD  transpiration  flux  density,  g  cm  " 2  s  " 1 

PT  potential  transpiration,  transpiration  calculated  with  the  max- 

imum input  k,,  cm3  h" 1 

R  total  water  flow  resistance  in  the  Soil-Plant-Atmosphere  Con- 

tinuum (SPAC),  s  cm-1 

DAYL  daylength,  seconds 

RH  relative  humidity,  % 
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APPENDIX  2.— Header  File  Format  and  Notes  for  H20TRANS  and  DAYTRANS. 

H20TRANS  Header  File 


Mnemonic 

Card  # 

Columns 

Format 

Description 

Default 

MAXX 

1 

1-2 

12 

NUMBER  OF  STATE 

- 

VARIABLES 

NEGX1 

3 

11 

OVERRIDE  NEC  STATE 

0 

VAR?  0  =  NO,  1  =  YES 

(DO  NOT  PRINT  DIAGNOSTIC) 

KSTART 

4-6 

13 

STARTING  TIME 

1 

KSTOP2 

7-10 

14 

ENDING  TIME 

365 

KSTEP 

11-14 

14 

TIME  INCREMENT(DT) 

1 

MAXB3 

16-17 

12 

HIGHEST  B  CONSTANT  IN- 

DEX 

- 

MAXG4 

19-20 

12 

TOTAL  #G  VALUES  TO  PRINT 

- 

MAXZ4 

22-23 

12 

TOTAL  #Z  VALUES  TO  PRINT 

- 

KPRINT 

25-26 

12 

PRINTER^  INTERVAL 

KSTEP 

KTHETA5 

28-29 

12 

PRINTER  STARTING  POINT- 

0 

INCLUSIVE  ITERATIONS 

FROM  BEGINNING 

NODUMP6 

1 

31-32 

12 

"l"  =  DO  NOT  DUMP 

0 

"0"  =  DUMP 

JCT 

2 

14 

A4 

CARD(DUMP)DECK  CODE 

"SOWR" 

T 

5 

9A8 

TITLE 

STEVE  RUNNING 

X(2,I]7 

3 

1-10 

E10.4 

INITIAL  CONDITION  STATE 

0.0 

VAR  #1  (8  PER  CARD) 

11-20 

E10.4 

INITIAL  CONDITION  STATE 

0.0 

VAR  #2  (8  PER  CARD) 

B(I)78 

(4) 

1-10 

E10.4 

B(l)  CONSTANT  (8  PER  CARD) 

0.0 

11-20 

E10.4 

B(2)  CONSTANT 

0.0 

IGP(I)9 

(5) 

1-2 

12 

INDEX  #  OF  G  VALUE  TO 

PRINT  (20  MAXIMUM) 

3-4 

12 

IZP(I)9 

(6) 

1-2 

12 

INDEX  3  OF  Z  VALUE  TO 

PRINT  (10  MAXIMUM) 

3-4  12 
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HEADER  FILE  NOTES 


1.  NEGX— If  this  code  is  "0,"  a  diagnostic  is  written  for  each  state  variable  that  becomes  negative  in  value.  To 
disable  this  diagnostic,  type  "1"  on  header  record. 

2.  KSTOP-For  hourly  resolution  model,  compute:  KSTOP  =  KSTART  +  (No.  days  ■  24). 

3.  MAXB— Always  enter  the  highest  index  number  of  B  constants  actually  used. 

4.  MAXG.  MAXZ— This  is  the  total  number  of  G  values  to  print. 

5.  KTHETA— This  feature  works  in  conjunction  with  KPRINT1.  Suppose  KSTART  =  200,  KSTOP  -  224,  KPRINT 
=  12,  and  KTHETA  =  0  (default).  Printing  will  occur  at  200,  211,  and  223  (i.e.  hours  1,  12,  24).  If  KTHETA  =  3, 
printing  occurs  at  202  and  214  (i.e.  hours  3, 15  and  24). 

6.  NODUMP— If  "0,"  unit  62  must  be  available  for  writing  upon;  if  "1,"  unit  62  is  not  used. 

7.  Card  #3  begins  initial  conditions  for  state  variables.  With  sequential  numbering,  1  through  8,  the  first  8  state 
variables  are  initialized.  On  card  #4,  #9  through  16  are  initialized,  etc.  Do  not  leave  fields  blank  until  the  MAXX 
(or  MAXB)  number  of  variables  has  been  reached.  Do  not  skip  any  fields,  as  the  processor  has  no  way  of  deter- 
mining the  index  number  of  the  initial  value  except  sequentially,  starting  with  1. 

8.  B  constants  are  handled  identically  as  X  initial  conditions.  Values  may  be  read  in  (punched)  as  floating  point  "F 
FORMAT."  E  FORMAT  is  not  mandatory.  However,  some  machines  require  E  FORMAT  numbers  to  be  punched 
in  a  specific  way.  Check  your  FORTRAN  reference  manual  at  your  installation  for  exact  details. 

9.  These  two  cards  contain  the  actual  index  numbers  of  the  variables  to  print  (G  and  Z  variables).  The  only  limita- 
tion is  that  they  be  less  than  or  equal  to  the  maximum  number  allowed  (20  and  10  respectively  for  G  and  Z).  The 
order  of  printing  is  as  follows: 

X(l)  THROUGH  X(MAXX)  -  STATE  VARIABLES, 

1(1)  THROUGH  I(MAXX)  -  FLOWS, 

G(i),  i  e  IG  =  (IGP(j),  j  =  1,  MAXG) 

Z(i),  i  e  Iz  =  (IZPQ),  J  -  1,  MAXZ) 
For  the  G's  and  Z's,  the  order  of  listing  is  the  order  of  printing,  (i.e.  there  is  no  resequencing  system  in  the  proc- 
essor to  reorder  the  G's  and  Z's  in  a  sequential  order  if  they  are  out  of  order.) 
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APPENDIX  3.— Input  Format  for  H20TRANS  Driving  Variables,  Z(i). 


Card         Columns        Format  Variable 


1-4 

A4 

Site  identification 

5-7 

A3 

Year 

8-10 

F3.0 

Julian  date 

12-13 

A2 

Card  ID 

14-17 

F3.1 

Daily  precipitation 

18-20 

F3.0 

Soil  temperature 

\A 

A4 

Site  ID 

5-7 

A3 

Year 

8-10 

A3 

Julian  date 

12-13 

A3 

Card  ID 

15-17 

F3.1 

Incoming  shortwave  radiation  (ly  min-1 

18-20 

F3.0 

Air  temperature  (°C) 

21-23 

F3.0 

Relative  humidity  (%) 

Card  2  format  is  repeated  on  cards  2  through  5  to  give  24  input  sets  per  day  on  a 
total  of  five  input  cards  per  day. 
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APPENDIX  4.— Input  Format  for  DAYTRANS  Driving  Variables,  Z(i). 


Card        Columns        Format  Variable 


\A 

A4 

Site  ID 

6-7 

A2 

Year 

9-11 

13 

Julian  date 

14-16 

F3.1 

Air  temperature  (°C) 

18-20 

F3.1 

Night  minimum  air  temperature  (°C) 

22-23 

F2.0 

Relative  humidity  (%) 

25-28 

F4.0 

Average  incoming  SW  radiation  (ly  min 

30-32 

F3.1 

Soil  temperature  (°C) 

34-37 

F4.2 

Precipitation  (cm) 

One  card  is  needed  for  each  day. 
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APPENDIX  5.— FORTRAN  Listing  of  H20TRANS. 

RE AL*8  F  N  A M 1 ,  F N  AM  2,  F N AM 3 
INTEGER  TTY 

COMMON  K,X(2,20),F(20,2C),G(100),B(100),Z(20),Y(160),MAXX,MAXB, 

1  KSTART,KSTOP,KSTEP,IGP(20), IZP(1 0),KPPINT,MAXG, MAXZ,T(20),JCT, 

2  NEGX, KTHETA, NOD U MP 

DATA  TTY/ IN, INI, I  OUT/ 5, 20, 21,  22/ 
WR 1TE( TTY, 5 ) 

5  FORM AT ( *  ENTER  HEADER  INPUT  FILE   NAME ( MAX   10  CHAR):  ',$) 
READ ( TTY, 6 )  FNAM1 

6  FORMAT ( Al 0 ) 
WRITE( TTY, 7) 

7  FORMAT ( *  ENTER  DATA   INPUT  FILE  NAMECMAX  10   CHAR):    ', $ ) 
READ (TTY, 6 )  FNAM2 

WRITE(TTY,6) 

8  FORMAT  (  *  ENTER  OUTPUT  FILE   NAME(MU  10   CHAR):  $) 
READ ( TT Y, 6 )  FNAW3 
0PEN(UNIT=IN1,FILE=FNAM1 ) 

OPEN (UNIT=IN, FILE =FNAM2) 
OPEN(UNIT  =  IO'JT,FILE  =  FNAM3) 
CPEN(UMT=23,FIL£  =  'PSN.DAT') 

C 

C  MAIN   PROGRAM,   CALLS  ALL  SUBROUTINES 

C  SHOULD  NOT  REQUIRE  CHANGE 

C 

CALL  HE AD ER( INI, TTY) 
K=KSTART 
10   CALL  ZCOMP(IN) 
CALL  PPOCES 
CALL  P.'tOTQ 

WRITE(23,100 )  Z(1),G(67) 
100  FORMATC1X,F5.0,1X,F6.1) 
CALL  FLOW 
CALL  YCDMP 
CALL  PRINT(IOUT) 
CALL  UPDATE 

IF   (K.GT.KSTOP-KSTEP)   GO  TO  20 
K?=K+KSTEP 
GO  TO  10 
20  K=KSTOP 
CALL  ZERO 
CALL  YCOMP 
CALL  PMNT(IOUT) 
STOP 
END 

SUBROUTINE  HEADER (INI , TTY) 
INTEGER  TTY 

COMMON  K, X( 2,20), F( 20,20 ),G (100 ),B( 100), Z( 20), Y( 160), MAXX, MAXB, 

1  KSTART,KSTOP,KSTEP, I GP (20 ) / I ZP ( 1 0 ) , KPR IN I , MAXG, MA XZ, T (2 0 ) , JCT/ 

2  NEGX,KTH  ET A, NOD U MP 
DIMENSION  IDUM(20) 

C 

C  READS  HEADER  INFORMATION  CARDS 

C  SHOULD  NOT  REQUIRE  CHANGE 

C 

READ (INI, 1000)   MAXX, NEGX, K ST ART, K STOP, K STEP, MA  XB, MAXG, MAXZ, 
1     KPRINT, KTHETA,NODUMP, JCT, T 
1000  F0RMAT(I2, II , 13, 2 14, 6 13,/, 2 1 A4) 

IF   (MAXX.LE.0.OR.MAXX.GT.20 )   CALL  ERROR (5, MAXX, TTY) 
IF   (MAXX. GT. 16)  GO  TO  20 
IF   (MAXX. GT. 8)   GO  TO  10 


30 


RE AD (IN  1/1001)   (X(2,I ),I=1,8) 

1001  FORMAT  (8E10.0) 
GO  TO  30 

10  READ(lNl,1002)   (X(2, I), 1=1, 16) 

1002  FORMAT   (8E10.0,/, 8E1 0.0) 

GO  TO  30 

20  READ(INl,1003)   (X ( 2, I), 1=1 , 20 ) 

1003  FORMAT   (  2 ( 8E 10 . 0, / ), 4E10. 0 ) 
30   IF   (MAXB.GT.32)  GO  TO  70 

IF   (MAXB.GT.24)  GO  TO  60 

IF   (MAXB.GT.16)   GO  TO  50 

IF   (MAXB.GT. 8)  GO  TO  40 

IF  (MAXB.EQ.0)  CALL  ERROR ( 1 , 0 , TTY ) 

READ(IN1,2000)   (B(I),  1=1,8) 

2000  FORMAT  (8E10.0) 
GO  TO  170 

40  READ(IN1,2001)   ( B ( I ), 1=1,1 6 ) 

2001  FORMAT   (8E10.0,/,8E1 0.0) 
GO  TO  170 

50  READ(IiNly  2002)   (B(I),  1=1,24) 

2002  FORMAT   ( 2 ( SE1 0 . 0, / ) , 8E10 . 0 ) 
GO  TO  170 

60  READ ( I Nl, 2003 )   ( B ( I ),  1=1,3 2 ) 
2C03  FORMAT  (3(8E10.0,/),8E10.0) 
GO  TO  170 
70  READ(IN1,2004)   ( B  ( I ),  1=1, 40 ) 
2004  FORMAT   ( 4( 8E10 .0, /), 8£10 .0 ) 
170   IF   (MAXG.EQ.O )  GO  TO  180 
J=21-MA.XG 

RE  AD  (I   1,3000)   (IGP(  I ),  1=1 ,  MAXG  ),  (  IDUM(  I),  1=1,  J) 

3000  FORMAT  (2112) 
180  J=11-MAXZ 

IF   (MAXZ.EU.O)   GO  TO  19C 

RE AD( INI, 3001 )  (IZP(I),I=1,MAXZ),(IDUM(I),I=1,J) 

3001  FORMAT  (1112) 

190   IF   (KSTART. EQ.U  )   KST ART=1 
IF  (KSTQP.EQ.O)  KSTOP=365 
IF   (KSTEP.EQ.O)  KST£P=1 
IF   (KPklNT.EQ.O)   KPRI NT=KSTEP 
IF   (JCT.EQ.*         ')  JCT='SCWR' 
DO  200  1=1,20 

IF   (T(l).NE.'         ')   GO  TO  210 
200  CONTINUE 

T(1)  =  *STEV 

T(2)='fc.  RU' 

T(3)='NNIN' 

T(4)=*G 
210  RETURN 

END 

SUBROUTINE  ZCOMP(IN) 

COMMON  K,XS(2, 20  )  , F ( 20, 20 ) , G ( 100 ) , B ( 1 00  ), Z  (  20  ) 

C 

C  REAUS  CLIMATIC  DATA  CARDS 

C  CHANGES   IN  CLIMATIC  DATA  INPUT  FORMAT   GO  HERE 

C  *     Z  FUNCTIONS  GO  HERE 

C  *     Z  DESCRIPTION 


C  -   

C  1  JULIAN  DAY 

C  2  PRECIPITATION  (CM) 

C  3  AIR  TEMPERATURE  (C) 
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C  4  RELATIVE   HUMIDITY   (PER  CENT) 

C  5  SOIL  TEMPERATURE  (C) 

C  6  RADIATION   (AVE.   LANGLEYS  PER  MIN) 

C  7  D  AYLENGTH   ( SUNRISE- SUNSET ) 

C  8  NIGHT  MINIMUM  AIR  TEMP.  (C) 

C 


100  F0RMAT(8X, I 3, F5. 1 , F4 . 1, F3. 0 , F5. 2 , F4. 1 , F5 .2 ) 

RE  AD  (IN/  100, END -50)   JD,  XT  AIR,  XTMI N,  XRHU*i,  XR  AD,  XTSO IL,  XPPT 

ZU)=JQ 

Z(2)=XPPT 

Z(3)=X?AIR 

ESD=6.1078*EXP( (17. 269*Z(3) )/ (237. 3+2(3) ) ) 

ES=XRHUM/100.*ESD 

VPD=AMAXI((ESD-ES),0.0) 
Z(4)=217.E-6*VPD/(Z(3)+273.16) 
Z(5)=XTS0IL 
Z(6)=XRAD 

XD=JD-79. 

IF(XD   • LT.   0.0)  XD=286.0+JD 

DAY=3.125*(SIN(XD*0.01721  ))+12.0 
Z(7)=DAY*3600*0.8 
Z(8)=XTMIN 
Z(12)=XRHUM 
Z(ll  )  =  VPD 
50  RETURN 
END 

SUBROUTINE  UPDATE (TT Y ) 
INTEGER  TTY 

COMMON  K, X(2,20),F(20,20),DUM(380 ), MAXX, DUMM( 47 ) , NEGX 

C 

C  UPDATES  STATE   VARIABLES    (COMPARTMENTS)   AFTER  EACH  ITERATION 

C  SHOULD   NOT  REQUIRE  CHANGE 

C 

DO  10   1=1, MAXX 
L=l 

X(2,I)=X(2,I)+X(1,I) 

IF   (X(2, I ) .LT.O.. AND. NEGX.EQ. 0)  CALL  ERROR ( 2,L, IT Y ) 
10  CONTINUE 
RETURN 
END 

SUBROUTINE  YCOMP 

COMMON  KyX(2,20),F(20,20),G(100  )y 3 ( 1 00 ) , Z ( 20 ) , Y ( 1 60 ) , M AXX, MAXB, 
1   KST AR T, KSTOPy  KSTEP, IGP ( 20 ) , IZP ( 1 0 ) , KPRINT, MAXG, MAXZ 

C 

C  PREPARES   MATRIX  OF   COMPUTED  VALUES  FOR  PRINTOUT 

C  SHOULD  NOT  REQUIRE  CHANGE 

C 

DO  10  1=1,MAXX 

J=I*MAXX 

Y(I)=X(2, I) 
10  Y(J)=X(1,I) 

MAX2=2*MAXX 

DO  20  I=1,MAXG 

J=MAX2*-I 

L=IGP(I) 
20  Y(J)=G(L) 

MAX3=MAX2+MAXG 

CO  30  i=l,MAXZ 

J=MAX3+I 

L=IZP(I) 
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30  Y(J)=Z(L) 
RETURN 
END 

SUBPOUTINE  ERROR ( I, J, TT Y  ) 
INTEGER  TTY 


C 

C  LIBRARY  OF  ERROR  STATEMENTS  CALLED  WHEN   MISTAKES  ARE 

C  MADE  IN  OTHER  PART  UF  THE  PROGRAM 

C  SHOULD   NOT  REQUIRE  CHANGE 

C 


IF   (I.LE.O   .OR.    I.GE.8)   GO  TO  99 
GO  TO  (1,2,3,4/5,6,7),! 

1  WRITE(TTY,1001 ) 

1001  FORMATC   ERROR   IN  B  CONSTANTS.  MAXB=0',/) 
STOP 

2  WRITE(TTY,1002.)  J 

1002  FORMATC   ERRHP  IN  SUb.  UPDATE  -  X ( I 2, ' ) • LT.O . *, / ) 
RETURN 

3  WRIT£(TTY, 1003)  J 

1003  FORMATC  ERROR  IN  SUB.  PROCESS.     GC,I2,*)    IS  LT  0.',/) 
RETURN 

4  WPITE( TTY,1004)  J 

10C4  FORMAT ( *  ERROR   IN   SPECIAL  FUNCTION   C,I2,  ")',/) 
RETURN 

5  WRITECTTY, 1005)  J 

1005  FORMATC   ERROR  IN  NO.   STATE  VARIABLES.  MAXX=",I3,/) 
STOP 

6  WRITE(TTY,1006) 

1006  FORMATC  ERROR   IN  YCOMP',/) 
RETURN 

7  WRITE(TTY,1007) 

1007  FORMAT ( *   ERROR  IN  SUB.  ZCPMP.',/) 
RETURN 

99  WRITE(?TY,9Q00)  I 
9000  FORMAT   (1H-, "ERROR  IN  SUB.   ERROR.     ERROR  CODE=',I3) 
END 

SUBPOUTINE  FLOW 

COMMON  K, XS( 2,20) ,F( 20,20 ),G( 100), B (100), Z(20),DUM( 160), MAXX 
DIMENSION  X(20) 

C 

C  COUPLES  FLOW  BETWEEN  COMPARTMENTS 

C  NEW  COMPARTMENTS  OR  NEW  COUPLING  BETWEEN  COMPARTMENTS   MUST  GO  HERE 

C 

DO  10  1=1,20 

DO  10  J  =  l,20 

F(I,J)=0.0 
10  CONTINUE 

DO  20   1=1, MAXX 

X(1)=XS(2,I) 
20  CONTINUE 

C 

C  *     FLO*   FUNCTIONS  GO  HERE 

C 

F<  1,  2)=G(50) 
F(  2,  2)=G(1) 
F(  2,  3)=G(53) 
F(  2,  5)=G<51) 
F(  2,  6)=G(52) 
F(  3,  4)=G(29) 
F(   7,  7)=G(9) 
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F<  8,  >j)=G(62) 
F(  9,  <})  =  G(66) 
F(10,10)=G(67) 
DO  30  I=1,MAXX 

XS(1,I)=0.0 

30  CONTINUE 

DO  50  I=1,MAXX 

DO  50  J=1,MAXX 

PART=F(  J,I)-F(I,J) 

IF   (I.EQ.J)  PART=F(I,I) 
50  XS(l,I)=XS(l,n«-PART 

RETURN 

END 

SUBROUTINE  PROCES 

COMMON  KyXS(2y 20 )yF (20,20 )/G( 100 ),B (100 )  ,  Z  (  20  )  ,  DUM(  1  60 ),MAXX 
DIMENSION  X(20) 

C 

C  AUXILIARY  EQUATIONS  USED   IN  THE  STRUCTURE  OF  THE  MODEL 

C 

IF   (  IFUG.NE.O)   GO  TO  20 

IFLAG=1 

DO  10  1=1,100 

G(I)=-1.0S32 
10  CONTINUE 
20  DO  30  I=1,MAXX 

X(I)=XS(2,I) 
30  CONTINUE 

C 

C  *     G  FUNCTIONS  GO  HERE 

C 

G(l  )=AMAX1((Z(2)-(B(4)/B(7))*B(6) )*3(7)/0.0) 

G(2)=X(2)/8(2) 

G(3)=X(3)/B(3) 

G(4)=AMAXl(B(8)*2(3)*a(7),0.0J 

1F(X(D    .LE.   0.0)  G(4J=0.0 
G(5)=B{4)/3(7) 
G(9)=Z(2)*B(7)-GU) 

G<10)  =  AMAXl(B(10>,0.2/G(2)«-0.01*rt(9),0.0) 

G(12)=10.*(1.0-EXP(-4.6*Z(6))) 

G(13)=1.-((G(5)-G(12D/G(5D 

IF(G(12)    .GT.   G(5))  G(13)=1.0 
G(14)=BC5)-CB(5)/(B(11)-B(10)))*(G(10)-B(10)) 

IF(G(10)   • GE.  3(11))  G(14)=0.005 
G(15)=G(14)-(G(14)*0.05*(Z(4)  *1 .OS+6-4. 0 ) ) 

IF(GClb)   .LE.   0.0)  G(15)=0.005 
G(16)=G(15)*G(15)*0.003*(Z(3)-10.0) 

1F(Z(8)   .LT.   0.0)  G(16)=AMAX1(G(15)+3.02*Z(8)/0.005) 
G(17)=G(16)*G(13) 
G(18)=Z(4)*G(17) 
G(19)=G(18)*B<4)*Z(7) 
G(20)=G(17)/B(5) 
G<21)=G<19)/B(7)*10. 
G(50)=G(4) 

IF(X(1)-G(4)   .LT.   0.0)  G(50)=X(1) 
G(51)=0.0 

IF((G(1)*G<4))/B(7)  .  GT.   BCD)     G(  5 1 )  =  (  G  (  50  )  -B  (  D  >  *  B  (  7  ) 
G(52)=0.0 

IF(X(2)*G(50)   .  GT.    3(2))     G ( 52 ) =X( 2 ) ♦ ;(5 0 )-G (51 ) -B ( 2 ) 
G(53)=o(3)-X(3)+G(19) 
G(55)=0.0 
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c 

C     PENMON  FUNCTION 
C 

G(28)=PENMON(Z(3),Z(6),Z(ll ),G(5),G(17) ) 
G(29)=G(28)*B(4)*Z(7) 

C 

RETURN 
END 

SUBROUTINE  PHOTO 

COMMON  K,X(2,20),F(20/20),G(100),B(100)/Z(20),Y(160) 

C 

C  Z(10)   =  CANOPY   AVERAGE  RADIATION 

C 

Z(10)=(Z(6)+Z(6)*EXP(-0.7*G(5)/2. ))/2. 

C 

C  Z(9)   =  AVERAGE   NIGHT  TEMPERATURE 

C 

Z(9)=(Z(3)   +  Z(B))/2. 

C 

G(64)=( (Z(IO)-O. 0143)/ (Z(10)+0. 322))* 
C(0. 0182+0. 0105*Z(3)-0. 000194* (Z( 3  )**2)) 

IF(G(64).LT.O. )G(64)=0. 

G ( 65 )  =  ( 0.000  6*(G(17)/1.6)*G(64) )/ 
C((G(17)/1.6)+G(64)) 

G(62)=0.001*(24.-Z(7)/3600.)*EXP(0.2*Z(9)) 

G(66)=G(65) *Z(7) 

G(67)=G(66)~G(62) 

RETURN 

END 

SUBROUTINE  ZERO 

COMMON  K,X( 2,20), F( 20,20 ), G (1 00 ), 8(100 ), Z( 20) 

C 

C  ZEROS  PARAMETERS    IN  MODEL 

C  SHOULD  NOT  REQUIRE  CHANGE 

C 

CO  20  1=1,20 

X(l, I)=0. 

Z( I)=0. 

DO  20  J=l,20 
20  F(i,J)=0.0 

DO  30  1=1,100 
30  G(I)=0.0 

RETURN 

END 

SUBROUTINE  PRINT(IOUT) 

COMMON  K,  K(  2,20),  F(  20/20  ),G  (100  ) ,  HQ 0 0  )  / Z (  20 ) /  Y  (  1 60  )  ,  M  AX X,  M AXB/ 

1  K.ST  ART,  KSTOP,KST£P,  I GP  ( 20  )  ,  1ZP  ( 1  0) /  KPP  INT ,  "A XG,  MA XZ,  T( 20  )  ,  JCT/ 

2  NEGX/KTHETA/ NODUMP 

DIMENSION  XNULL(100),LBL(100,2), JK(10) 
DATA  JK./ 1,8/  15,22,29,36,43/50,57/  64/ 


C 

C  LABELS  LINE  PRINTER  OUTPUT 

C  ESTABLISHES  OUTPUT  FORK AT 

C  PRINTS   OUT   ALL   REQUESTED  PARAMETERS 

C  SHOULD  NOT  REQUIRE  CHANGE 

C 


IF   (IFLAG.EQ.l)   GO  TO  30 

1FL AG=1 

DO  10  1=1,100 

XNULL(1)=0.0 
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DO  10  J=l,2 
10  LBL(I,J)='  * 
J=100-MAXB 

WRITE  (IOUT/1000)   T, M AXX, ( B ( I ), I=1,MAXB ), 
1  (XNULL(I),i=l,J),XSTART,FSTOP,KSTEP 
1000  FORMATC '1 *,/, '1 ',/,5X,'KXBZ   SIMULATION  ? ROC ESSDR * , /, 5X/20 A4, 

1  5 X, 'STATE  VARI ABLES=', 12, /,5X, 'B  CC)NSTANTS  =  ',/, 

2  10(5X,10(1X,1PG11.4),/),5X,  'INITIAL  THE ',  3X,  13, /,5X, 

3  'TE RMINATIQN',4X, 13,/, 5 X, 'INCREMENT  (DT)',1X,I3) 
MAX2=2*MAXX 

MAX3-MAX2+MAXG 

MY=MAX3+MAXZ 

JJ=100-MY 

IF  (MY.LE.50)  JJ=JJ-50 
DO  20  1=1, MY 

IF   (I.LE.MAXX)  L3L(I/1)='X( ' 
IF   (  I.GT.MAXX. AND.I.LE.MAX2)   LBL( 1,1 )=*I{ ' 
IF   (I.GT.MAX2.AND.I.LE.?4AX3)   LBL(  1,1)  =  ' J (  * 
IF   (1.GT.MAX3.  AND. I. LE. MY)   LBL(  I,  1  )■=  'Z(  * 

20  CONTINUE 

DO  21   1=1, MUX 
ENCODE(3,2002,LBL(I,2))  I 
20  02  FORMAT { 12,')') 

21  CONTINUE 

DO  22  I=1,MAXX 
J=MAXX+I 

£NCODE(3,2002,LBL(J,2))  I 

22  CONTINUE 

DO  23  I=1,MAXG 

J=MAX2+I 

L=IGP<I) 

ENCODE(3,20Q2,LBL(J,2 ))  L 

23  CONTINUE 

DO  24  I=1,MAXZ 

J=MAX3+I 

L=IZP(I) 

ENCODE(3,2002,LBL(J,2))  L 

24  CONTINUE 

30  IF   (KFLAG.EQ.l)  GO  TO  40 
KFLAG=1 
ICT  =  0 

WRITECIOUT/2000)  T 

2000  FORMAT ( *1',/,5X,20A4,  /  ) 
WRITE(IOUT,2001 )    ((LBL{ I, J), J =1,2), 1=1, 100) 

2001  FORM  AT ( '         TIME ', 10 ( 7X, A2, A3 ), /, 9 ( 8 X, 1 0( 7 X, A2, A3 ) // ) ) 
ICT=ICT+15 

40   IF   ( MOD (K-KST ART  +  1  +  KPRINT-K  THETA, KPRINT)«NE.O)   GO  TO  60 
IF   (MY.LE.50)  GO  TO  50 

WRITE  (IOUT/3000)   K, ( Y ( I )/ I =1, MY ), ( XN ULL ( I ) ,1 =1, JJ) 
300 J  FORMAT(3X,I4,10(1X,1PG11.4),//9(7X,10(1X,1P511*4),/)) 
ICT=ICT+11 
GO  TO  60 

50  WRITE  CIOUT/3001)   K, ( Y ( I ),  I =1 , MY ) , ( XNULL ( I ) , I =1 , J J ) 
3001  FORMAT<3X,I4,10ax,lPGll.4),/,4<7X,10ClX,lPG11.4),/)) 

ICT=ICT+6 
60  IF   (ICT.GT.50)  KFLAG=0 

IF   (NODUMP.SQ.l )  RETURN 

UPhN(U*IT=2  3,FILE='DUMP.DAT') 

KK=0 

DO  70  1=1, MY/7 


36 


J1=J 

J2=I+6 
KK=KK+1 

WRITE (23/ 40  0  0 )  JCT/K^JK(KK)/ (Y(L),L=J1/J2) 
4000  FORMAT  ( A4,  I4,1X, I2/7E10.3) 
70  CONTINUE 
RETURN 
END 

FUNCTION  PENMON(T AIR, RAD, VPD, XL  A  I , CQND ) 

GAMVA=0. 646*0. 000 6*T AIR 

PLAI=XLAI/2. 

T1=TAIR>0.5 

T2-TAIR-0.5 

SVPl=6.1078*£XP((17.269*Tl)/(237.0«-Tl)) 

SVP2=6.1078*EXP<(17.269*T2)/(237.0+T2J ) 

SL0PE=SVP1-SVP2 

XNETR=RAD*0. 8*697.3 

CP=l-0lE+3 

PA=1. 292-0.  00428*TAIP. 
RA=5.0 

RS=100./COND 

XLAT=(2.501   -  U.U024*TAlk)   *  1.02+6 
XTRANS=( (SLJPE*XNETR)/PLAI*(CP*PA)*( VPD/RA) )/ 
C(SLOPE*GAMMA*(1.0+RS/RA)) 
PENMON=XTRANS/(XLAT  *  10.) 
RETURN 
END 

FUNCTION  S4(K) 

C 

C  SETS  MODEL  STARTING  POINT  TO   REQUESTED  JULIAN  DATE 

C  REGARDLESS  OF  DATA. 

C  SHOULD  NOT  REQUIRE  CHANGE 

C 

IN  =  20 

10  READ(IN/1000/END=20 )  KD 
100  0  FORMAT (  7X,I3) 

IF   (K.GT.KD)   GO  TO  10 
BACKSPACE  IN 
S4=l. 
RETURN 
20   REWIND  IN 

IFLAG=IFLAG«-1 

IF   (IFLAG.Ew.2)  STOP 

GO  TO  10 

END 
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C  DATA  FILE   INPUT  STATEMENTS 

C 

RE AL*8  F  N  AMI , F N AM 2,  F  N  AM3 
INTEGER  TTY 

COMMON  K,X(2,20),F<20,2C),G(100)/3(100),Z(20),Y(160),MAXX,MAXB, 

1  KSTAPT, KST9P,KSTEP,IGP(20),IZP(10),KPRINT,MAXG,MAXZ,T(20),JCT, 

2  N  EGX, KTHET  A, NODUMP 

DATA  TTY/ IN, INI, I  OUT/ 5, 20, 21, 22/ 
WRITE(TTY,  b ) 

5  FORMATC  ENTER  HEADER  INPUT  FILE  NAME  ( MAX  10  CHAR):  ',$) 
READ  CTTY, 6 )  FNAM1 

6  FORMAT(AIO) 
WRITE( TTY,7) 

7  FORMAT ( *  ENTER  DATA   INPUT  FILE  N  A  M  E ( M  A  X  10   CHAR):    *,  $  ) 
READ ( TTYy  6 )  FNAM2 

WRITE(TTY,8) 

8  FURM AT ( *  ENTER  OUTPUT  FILE   NAME (MAX  10  CHAR):        $ ) 
READ ( TTY, 6 )   FN AM 3 

UPEN(UMT=IN1,FILE  =  FNAM1) 
OPEN  (  UN  IT  =  IN, FILE -=FNAM2) 
OPEN (UN  I T= I3UT, FI LE  =  F  NAM3 ) 

C 

C  MAIN  PROGRAM,   CALLS  SUBROUTINES 

C 

CALL  HEAD£R(IN1,TTY) 
K  =  KSTArlT 
10   CALL  ZCOMP 
CALL  PEOCES 
CALL  FLOW 
CALL  YCOMP 

IF(Z(6).EQ.   0.0)   GO  TO  15 

CALL  PRINT(IOUT) 
15  CONTINUE 

CALL  UPDATE(TTY) 

IF   (K.GT.KSTOP-KSTEP)   GO  TU  20 

K  =  K>KSTEP 

GO  TO  10 
20  K=KSTOP 

CALL  ZERO 

CALL  YCOMP 

CALL  PHINT(IQUT) 

STOP 

END 

C 

SUB ROUT  IN  E  HEADER (INI, TTY) 

COMMON  K, K( 2,20), F( 20,20 ),G (100 ),B( 100), Z( 20), Y( 160), MAXX,MAX8, 

1  KST APT, KSTOP,KSTEP,IGP(2C),IZP(10),KPRINT,VAXG, MAXZ,T(20),JCT, 

2  NEGX,KTHET A, NODUMP 
DIMENSION  IDUM(20) 

READ  ( IN  1,  1000  )   MA XX,  NEGX, KST ART, K STOP, <ST EP, MA XB, M AXG, MA XZ, 
1  KPRINT,KTHETA,NODUMP,JCT,T 
10C0  FORMAT( 12, II, 13, 214, 613,/, 21 A4) 

IF   (MAXX.LE.O   .OR.   MAXX.GT.20)   CALL  ERRDR  (5,MAXX,TTY) 
IF   (MAXX.GT.16)   GO  TO  20 
IF   (MAXX.GT.8)   GO  TG  10 
READ(IN1,1001)   (X(2, I), 1=1,8) 

1001  FORMAT  (8E10.0) 
GO  TO  30 

10  READ(IM,1002)  (X(2,  I), 1=1, 16) 

1002  FORMAT   (8E 10.0,/, 8E1 0.0) 
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GO  TO  30 

20  READ(IN1,1003)   (X(2, I ), 1=1,20  ) 
10C3  FORMAT  (2(8E10.0,/),4E10.0) 
30   IF   (MAXB.GT.32)  GO  TO  70 
IF   (MAXB.GT.24)  GO  TO  60 
IF   (MA  KB. GT. 16 )  GO  TO  50 
IF  (MAXB.GT.8)  GO  TO  40 
IF  (MAXB.EQ.03  CALL  ERROR(1,0,TTY) 
READ(IN1,2000)   (B(I), 1=1,8) 

2000  FORMAT  (8E10.0) 
GO  TO  170 

40  READ(IN1,2001)   (  6  {  I ),  1=1,1 6) 

2001  FORMAT   (8E10.0,/, 8E10.0) 
GO  TO  170 

50  READ(I>ll,  2002)   (B(I),  1=1,24) 

2002  FORMAT   (2( 8E10. 0,/ ), 8E10.0) 
GO  TO  170 

60  READ(IN1,2003)   (8 (I ), 1=1,32) 

2003  FORMAT  (3( oSlO.O, /),8E10.0) 
GO  TO  170 

70  READCIN1/2004)   (B (I ), 1=1, 40 ) 

2004  FORMAT   ( 4( 8E10.0, /), 8E10.0 ) 
170  IF  (MAXG.EQ.O)  GO  TO  180 

j=21-MAXG 

READ ( I Nly 30  00)   (  I GP  (  I  ),  1=1 ,  M  AXG  ), ( IDUM( I), 1=1, J) 

3000  FORMAT  (2112) 
180  J=11-MAXZ 

IF  (MAXZ.EQ.O)  GO  TO  190 

READ( 1^1,3001)   (IZP(I),I=1,MAXZ), (IDUM( I),I=1,J) 

3001  FORMAT  (1112) 

190   IF   ( KSTART. EQ.O )   KST ART=1 
IF   (KST'IP.EQ.O)   KST0P  =  365 
IF   (KSTSP.E^.O)  KSTEP=1 
IF   (KPKINT. EQ.O )   KPR I NT=KSTEP 
IF   (JCT.EQ.4H         )  JCT=4KSOWR 
DO  200  1=1,9 

IF  (T(l).E(J.BH  )   GO  TO  200 

GO  TO  210 
200  CONTINUE 

7(1)=*STEV 

T(2)='E  RU* 

T(3)='NNIN' 

T(4)="G 
210  CONTINUE 

RETURN 

END 

C 

SUBROUTINE  ZCOMP 

COMMON  K,XS(2,20),F(20,20),G(100),B(100),Z(20) 


c 
c 
c 

*  z 

DESCRIPTION 

c 

1 

JULIAN  DAY 

c 

2 

PRECIPITATION  (CM) 

c 

3 

AIR  TEMPERATURE  (C) 

c 

4 

RELATIVE  HUMIDITY  (%) 

c 

5 

SOIL  TEMP  (C) 

c 

6 

RADIATION   (LANGLEYS   PER  A I M ) 

c 

7 

HOUR  OF  DAY 

DIMENSION  VAL(10,24) 
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IF   (IFLAG.NE.O)  GO  TO  10 

IFLAG=1 

XI=S4(K) 

IF   (XI.NE.l.)  STOP 
10  CONTINUE 

C 

C  *     Z  FUNCTIONS  GO  HERE 

C 

XI  =  S2(  VAL  (1,1)  ) 

Z(l)  =  VALU/l) 

Z(2)  =  VAL(2,1) 

Z(3)  =  VAL(5,1) 

Z(4)   =  VALC6/1) 

Z(5)  =  VAL(3,1) 

Z(6)  =  VAL(4/1) 

ZC7)   =  Z(7)   ♦  1.0 

IF   (Z(7)   .GE.   25.0)     Z(7)  =  1.0 

RETURN 

END 

C 

SUBROUTINE  UPDATE (TTY ) 
INTEGER  TTY 

COMMON  K/  X(2/20)/F(20/20)/DUM(3  80)/MAXX,DUMM(47),NEGX 

DO  10  I=1/MAXX 

L=I 

X{2/I)=X(2/I)+X{1,I) 

IF   (X(2,I) .LT.O.    .AND.   NEGX.EQ.O)   CALL   ;:RROR ( 2,L, TT Y ) 
10  CONTINUE 
RETURN 
END 

C 

SUBROUTINE  YCOMP 

COMMON  K, X(2,20),F(20,20),G(10Q ), 3(1 0  0 )/ Z ( 20) , Y ( 160 ), MAXX, MAXB, 
1  KSTART,KSTOP,KSTEP,IGP(20),IZPU0)/KPRINT/MAXG/MAXZ 

DO  10  1=1/MAXX 

J=I+MAXX 

Y(I)=X(2,I) 

Y(J)-X(1/ I) 
10  CONTINUE 

DO  20  l=l,M\XG 

J=2*MAXX+I 

L=IGP(l) 

Y(J)=G(L) 
20  CUNTINUE 

DO  30  I=1,MAXZ 

J=2*MAXX+MAXG+I 

L=IZP<I) 

Y(J)=Z(L) 
30  CONTINUE 

RETURN 

END 

C 

SUBROUTINE  ERROR   (I, J, TTY) 
INTEGER  TTY 

IF   (I.LE.O    .OR.    I.GE.8)   GO   TO  99 
GO  TO  (1/2/3/4,5/6,7)/! 

1  WPITE(TTY,1001) 

1001  FORM AT ( 1  HO / *  ERROR  IN  3  CONSTANTS.  MAXS=0',/) 
STOP 

2  1*RITE(TTY,  1002)  J 


40 


1002  FORMAT   (1H0, 'ERROR 
RETURN 

3  WRIT£(TTY/1003)  J 

1003  FORMAT  (1HO,'ERROR 
RETURN 

4  WRITE(TTY,1004)  J 

1004  FORMAT    (1HO, 'ERROR 
RETURN 

5  URITE(TTY,1005)  J 
ICO1)  FORMAT   (1  HO  /  'ERROR 

STOP 

6  WPITE(TTY,1006j 
100b  FORMAT   (1H0, 'ERROR 

RETURN 

7  KRITE(TTY,1007) 
1007  FORMAT   ( 1  HO /  'ERROR 

RETURN 
99  WPITE(TTY,9000)  I 
9000   FORMAT   (1H-, 'ERROR 
END 

C 

SUBROUTINE  FLOW 

COMMON  K, XS(2,20),F(2C, 20 ), G ( 100 ), B( 1 0  0 ), Z ( 20 ) ,DUM( 1 6  0 ), MAXX 

DIMENSION  X(20) 

DO  10  1=1,20 

DO  10  J  =  l,20 

F(I,J)=0.0 
10  CONTINUE 

DO  20   1=1, MAXX 

X(1)=X3(2,I) 
20  CONTINUE 

C 

C  *     FLOW  FUNCTIONS  GO  HERE 

C 

F(3,3)=G(50)-G(53) 

F(1,3)=G{51) 

F(3,2)=G(54) 

* (3,4)=G(55) 

F(4,5)=G(56) 

F(5,6)=G(57) 

F(d,9)=G(18) 

F(7,8)=G(62) 

F(4,7)=G(67) 

F(5,7)=G(68) 

F(10/10)=G(53)  ♦  G(58) 

DO  30  1=1, MAXX 

XS(1,I)=0.0 
30  CONTINUE 

DO  50   1=1, MAXX 

DO  40  J=1,MAXX 

IF   (I.EQ.J)   GO  TO  60 

XS(1,I)=XS(1,I)*F(J/I)-F(I/J) 
40  CONTINUE 
50  CONTINUE 

RETURN 

bO   XS(l,I)=XS(l,i;*F(I,I  ) 

GO  TO  40 
END 

C 

SUBROUTINE  PROCES 


IN  SUB.    UPDATE  -  X('/I2/').LT.O.',/) 

IN  SUB.   PROCES.     G(',I2,')    IS   LT   0 . ', / ) 

IN  SPECIAL  FUNCTION  (',12,')',/) 

IN  NO.  STATE  VARIABLES.  MAXX='/I3,/) 

IN  YCOMP',/) 

IN  SUB.   2C0MP. ',/) 

IN  SUB.   ERROR.     ERROR  C0DE=',I3) 
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COMMON  K,XS(2,20)/F(20,20)/G(100),b(100)/Z(20),DUM(160),MAXX 

DIMENSION  X  C  20  J 

IF   (IFUG.NS.O)   GO  TO  20 

IFLAG=1 

DO  10  1=1,100 

G(1)=-1.0E32 
10  CONTINUE 
20  DC  30  I=1,MAXX 

X(I )=X5(2y  I) 
30  CONTINUE 

C 

C  *     G  FUNCTIONS  GO  HERE 

C 

G(1)=SJ(Z(3),Z(4)) 

G(2)=B(9)/B(10) 

G(3)=X(3)/6(3) 

G(4)=X(4)/3(4) 

G<5)=X(5)/B(5) 

G(6)=(C(4)+G(5))/2. 

G(7)=X(7)/3(7) 

G(10)  =  AMAXl(B(13)y0.2/G{6)*0.01*8Ql),0.0) 
G(11)=10.*(1.0-EXP(-4.6*Z<6))) 
G(12)=1.0-(G(2)-G(11))/G(2) 
IF(G(11).GT.G(2))  G(12)=1.0 

GU3)  =  AMAX1(3(12)-(BU2)/(B(14)-3(13)) ) * ( G( 10 )-B ( 13) ) , 0. 00 5) 

IF(G(10).LT.B(13))  G(13)=E(12) 

G(14)=G(13)-(G(13)*.05*(G(l)*1.0E6-4. ) ) 

IF(G(14)   .LE.   0.0)  G(14)=0.005 

G(15)=S1(Z(3)/Z{7)/G(14)) 

G(16)=G(15)*G(12) 

G(17)=G(16)*GCD 
G(13)=G(17) *B(9)*3600. 

C 

IF(Z(7)   .GT.  4.)Z(2)=0.0 

G<50)=AMAXl((Z(2)/4.-B(6)*(B(9)/y(10)))*B(10),0.0) 

G(51  )=AMIN1(Z(3)*9(1)*B(10)/X(1 )) 

IF(Z{3).L£.0.0)  G(51)=0.0 

G(52)=EXP(5.*(G(3)-1.0)) 

G(53)=G(52)*G(1)*3(10 )*360  0. 

IF(G(3).LE.0.O)  G(53)=0. 

IF(X(l).GT.O.)  G(53)=0. 

G{54)=AMAX1((G(50)+G(51)-(B(2)*B(10))),0. ) 

G(55)=AMAXl((G(5O)*G(51)-G(53)-G(54))-(B(3)-X(3)),0.0) 

G(56)=AMAX1(G<55)-<B(4)-X(4)),0.0) 

G(57)=AMAX1 (G(56)-(8(5)-X(5)),0.0) 

G(58)  =  U(2)/4.)*B(10)-G(50) 

C 

G(60)  =  B(8)*(1.-GU0)/B(14)  +  B(13)/B(14)) 

G(61)=X(8)/B(16) 

G(62)=U. 

IF<G(61).GE.G(18))   G( 62)-AM AXI (G ( 60)-X(3)/G(18)/B(16)) 
IF<G(61).LT.G{18))  G ( 62)=G ( 1 8)-G( 61 ) 
G(66)=X(7)/B(17) 

G(b7J)  =  AMINl(G(62)-G(b6),<B(7)-X(7))/B(17)/X(4)/B(17)) 
G(68)=AMIN1(G(62J-G(66)-G(67)/ ( B ( 7 ) -X ( 7 ) ) / 8 ( 17 ) ) 
G(69)=G(67)«-G(68) 

C 

IF(Z(5)   .LT.   0.0)G(70 )=10.Q 
IF(Z(5)   .GT.  0.0)  G{70)=1.0/Z(5) 
G<71  )=-0. 28+2.*G(10) 
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G(72)  =  G(70)«-G(71) 

G(73)  =  AMINKG(10)+G(72)*(G(69)*278./B(9)),B(14)) 

G<76)=G(lb)/B(12) 

G(77,=G(18)/B(10)*10. 

G(78)=0.0 

RETURN 

END 

SUBROUTINE  ZERO 

COMMON  K/X(2/20)/FC20/20)/G(100)/B(100)/Z(20) 
DO  10  1=1,20 

X(ly I)=0.0 

Z(I)=0.0 
10  CONTINUE 

DO  20  1=1/20 

DC  20  J  =  l,20 

F(I,J)=0.0 
20  CONTINUE 

DO  30  1=1,100 

G(I)=0.0 
30  CONTINUE 

RETURN 

END 

SUBROUTINE  PRINT(IOUT) 

COMMON  K, X( 2,20)/ F( 20,20 ),G(100 ), B(100) ,Z( 20), Y( 160) ,MAXX,MAXB, 

1  KSTART,KSTOP,KSTEP, I GP { 20 )  ,  IZP ( 1 0 ), KPR INT, ^AXG, M A XZ, T ( 2 0) , JCT, 

2  NEGX,KTHETA,NODUMP 

DIMENSION  XNULL(IOO), LBL(10  0,2), JK(10) 
DATA  JK/ 1,8, 15,22,29,36,43,50,57,64/ 
IF   (IFLAG.EQ.l)  GO  TO  30 
IFL AG=1 
DO  10  1=1,100 
XNULL(I)=0.0 
DO  10  J=l,2 
10  LBL(I,J)='  ' 
J-100-MAXB 

WRlT£(lOUT,10  00)   T,MAXX, (B ( I ), 1=1 , MAXB ) , 
1  (XNULL(I),I=l,J),KSTART,KSTOP,KSTEP 
1000  FORMAT  ( '1 ', /, '1 ', /, 5  X, 'KXBZ   SIMULATION  P ROCESSOR ',/, 5X, 20 A4, /, 

1  5X, 'STATE  VARIA8LE5= ', 12, /, 5X, 'B  CONST ANTS= ', /, 

2  10(5X,10(1X,1PG11. 4), /),5X, 'INITIAL   TI ME' , 3X, I  4, /, 5 X, 

3  'TERMINATION', 4X, 14, /,5X, 'INCREMENT  (0T)',1X,I4) 
MAX2=2*MAXX 

MAX3=MAX2+MAXG 

MY=MAX3+MAXZ 

JJ=100-MY 

IF   (MY.LE.50)  JJ=JJ-50 
DO  20  1=1, MY 

IF   (I.LE.MAXX)   LBL( I,1)  =  *X(  ' 
IF   (I.GT.MAXX.AND.I.LE.MAX2)   LBL (  I, 1 )  =  '  I< ' 
IF   (I.GT.MAX2. AND.I.LE.MAX3)   LBL (  1, 1 )  =  ' G( ' 
IF  ( I.GT.MAX3. AND. I. LE. MY)   LBL( I, 1  )= 'Z( ' 

20  CONTINUE 

DC  21  I=1,MAXX 
ENCODE   (3,2002,LBL(I,2))  I 
2002  FORMAT (12, ' ) ' ) 

21  CONTINUE 

DO  22  I=1,MAXX 
J=MAXX+I 
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ENCODE  (3,2002,LBL(J,2))  I 

22  CONTINUE 

DO  23  I=1,MAXG 

J=MAX2+I 

L=IGP(I) 

ENCODE  (3,2002,LBL(J,2))  L 

23  CONTINUE 

DO  24  I=1,MAXZ 

J=MAX3+I 

L=IZP(I) 

ENCODE  (3,2002,LBL(J,2))  L 

24  CONTINUE 

30  IF  {KFLAG. EQ.l )  GO  TO  40 
KFLAG=1 
ICT=0 

fcRIT£(lOUT,20OO)  T 

2000  FORMAT(*1',5X,20A4,/) 
WRITE(IOUT,2001)   (UBLU/J),  J  =  l,2),  1  =  1,  100) 

2001  FORMAT   (4X, 'TIME ', 10 ( 3X, A2, A3, 4X) ,// 9( 3X, 1 0 ( 3X, A2, A3, 4X) ,/ ) ) 
ICT=ICT*13 

40  IF   ( MOD(K-KST ART+l+KPRINT-KTHET  A, KPRINT).NE.O)  GO  TO  60 
IF   (MY.LE.50)  GO  TO  50 

WRITEC  iOUT/3000)  K, ( Y ( I ), 1  =  1 ,MY), ( XNULLC I) , 1=1 , JJ ) 

3000  FORMAT  (3X,I4,10(1X,1PG11.4),/,9(7X,10(1X,1PG11.4),/)) 
ICT=ICT*11 

GO  TO  oO 

50  WRITE(IOUT/3001)   K, ( Y ( I), I =1, MY ) / ( XNULL( I) , I =1/ JJ ) 

3001  F0RMAT(3X, I4,10(1X,1PG11.4),/,4(7X,10(1X,1PG11.4),/)) 
ICT=ICT*6 

oO  IF  (ICT.GT.50)  KFLAG=0 
IF   (NGDUMP.SQ.l)  RETURN 
GPEN(UNIT=2  3,FIL£='DUMP.DAT*) 
KK=0 

DO  70  1=1, MY, 7 
J1  =  I 
J2=I+6 
KK=KK+1 

WRITE  (23,4000)  JCT,K, JK(KK), ( Y ( L), L=J1, J2) 
4000  FORMAT  {  A4,  1 4,  1  X,  1 2,  7 E10 .  3  ) 
70  CONTINUE 
RETURN 

END 

C 

FUNCTION  51 (rAIR,XHOUR/COND) 
IFCXHOUR-5.)  91,90,91 

90  SRT=TAIR 

91  IF(SRT.GT.O)   GO  TO  92 
Sl=AMAXl(CaND+O.O2*SRT/0.0  0  5) 
IFCSRT.LT.-7.)  Sl=0.005 

GO  TO  93 

92  Sl=CONu+cnND*0.00  3*(TAIR-10.0 ) 

93  RETURN 
END 

C 

FUNCTION  S2  (VAL) 
DIMENSION  VAL(10,24) 
IN  =  20 

IF   (IFLAG.EQ.l)   GO  TO  20 

READ  (IN, 1000)  CVAL(I,1),I-1/3),C(YAL(I,J),I=4,6J,J=1,6), 
1  ((VALCI/J)/I=4,b),J=7,12)/((VAL(I/J),I=4,6),J=13,18), 
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2  ((VAL(I,J),I=4,6),J=19y24) 
1000  FORMAT  ( 7X, F3 .0, 4 X, F3 . 1, F3- 0/ /, 3 ( 1 4X, b( F3 . 1 , 2F3 . 0 ) , / ), 
1  14X,6(F3.1,2F3.0 )) 
6  K=l 

IFLAG=1 

DO  10   1  =  2,24 

VAL(1,I)=VAL(1,1J 

VAL(2,I)=VAL(2,1) 

VAL(3,I)=VAL(3,1) 
10  COM  INt/E 

S2=l<. 

K  =  K  +  1 

RETURN 
20  DO  40  J=l,10 

DO  30  1=2,24 

VAL(J,1-1)=VAL(J,I) 
30  CONTINUE 

VAL(J,24)=0.0 
40  CONTINUE 

IF  (K.EQ.24)   IFLA G=0 

K  =  K*1 

S2=2. 

RETURN 

END 

C 

FUNCTION  53(TAvRH) 
C  ******ABSOLUTE  HUMIDITY  DFFICIT   (AIR  TEMP,    REL  HUM) 

ESD  =  6.1078  *  EXP  ((17.269  *  Ik) /  (237.3  +  T A ) ) 
ES=RH/100. *E3D 

c  ******       £N  XHE  EVENT  DE5*   IS  GTP   THAN   MR,    RETURN  ZERO 

VPD  =  AMAX1 ( (ESD  -  ES),  0.0) 
53=  217. E-6  *  VPD  /   (TA  +  273.16) 
RETURN 
END 

C 

FUNCTI'JN  S4(K) 
IN  =  20 

10  READ   (IN.,1000,END=20  )  KD 
1000  FORMAT (  7X,I3) 

15  IF   (K.GT.KD)  GO  TO  10 

BACKSPACE  IN 

S4=l. 

RETURN 
20  REWIND  IN 

IFLAG=1FLAG+1 

IF   (IFLAG.EQ.2)  STOP 

GO  TO  10 

END 
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U.S.  Department  of  Agriculture 
Forest  Service 

Rocky  Mountain  Forest  and 
Range  Experiment  Station 

The  Rocky  Mountain  Station  is  one  of  eight 
regional  experiment  stations,  plus  the  Forest 
Products  Laboratory  and  the  Washington  Office 
Staff,  that  make  up  the  Forest  Service  research 
organization. 

RESEARCH  FOCUS 

Research  programs  at  the  Rocky  Mountain 
Station  are  coordinated  with  area  universities  and 
with  other  institutions.  Many  studies  are 
conducted  on  a  cooperative  basis  to  accelerate 
solutions  to  problems  involving  range,  water, 
wildlife  and  fish  habitat,  human  and  community 
development,  timber,  recreation,  protection,  and 
multiresource  evaluation. 


RESEARCH  LOCATIONS 

Research  Work  Units  of  the  Rocky  Mountain 
Station  are  operated  in  cooperation  with 
universities  in  the  following  cities: 


Albuquerque,  New  Mexico 

Flagstaff,  Arizona 

Fort  Collins,  Colorado* 

Laramie,  Wyoming 

Lincoln,  Nebraska 

Rapid  City,  South  Dakota 

Tempe,  Arizona 


'Station  Headquarters:  240  W.  Prospect  St.,  Fort  Collins,  CO  80526 


